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 ‡ 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 §. 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 at horizon crossing during inflation, the primordial bispectrum is dominated by equilateral configurations (k 1 k 2 k 3 ). Non-canonical kinetic terms will also yield non-Gaussianities of this shape.
The non-Gaussianities produced by the most general single-field inflation models with approximate shift symmetry have shapes that vary from being peaked on equilateral configurations to being peaked on enfolded (k 1 k 2 k 3 /2) configurations [13]. These types of non-Gaussianities, have been shown [13] to be generically well described by a linear combination of two bispectrum shapes, one where the bispectrum signal peaks on equilateral configurations and another one, orthogonal to it, called orthogonal shape.
Finally, modifications of the initial-state of the inflaton field, leave their signatures in a bispectrum dominated by flattened or enfolded configurations [7,14].
The standard observables to constrain primordial non-Gaussianity are the cosmic microwave background (CMB) and large-scale structure. Both CMB and largescale structure can measure, in principle, the bispectrum shape dependence and thus discriminate the shape of non-Gaussianity. The large-scale structure bispectrum however, includes a large contribution form non-linear gravitational evolution and biasing; compared with this contribution, the primordial signal "redshifts away" during the Universe's evolution [15,16,17]. On the other hand large-scale structure can probe scales non easily accessible from the CMB (e.g., see discussion in [18,19]) and offers other powerful probes. One technique is based on the abundance of rare events as they trace the tails of the underlying distribution [20] and has received renewed attention e.g., [18,21,22,23,24,25,26] and references therein). Recently, the effect of non-Gaussian halo bias [27,28] has been demonstrated to be an extremely promising avenue ‡ There are some cases where the trispectrum may be important (when, for example, the bispectrum is zero) but, in general, one expects the trispectrum contribution to be sub-dominant compared to the bispectrum one. § This is also the non-Gaussianity of standard slow-roll inflation, but, as mentioned above, in this case the amplitude is unmeasurably small, but see e.g., [12].
both for present [29,30,31] and future [32,33] data. This approach uses the fact that the power spectrum of density extrema (where galaxies are formed) on large scales increases (decreases) for a positive (negative) f NL . The signal, for these two techniques is given by an integral over the bispectrum. The shape-dependence is thus indirect and suppressed, but signatures still remain [23,12].
While techniques have been explored and developed to create non-Gaussian CMB maps with given bispectrum and power spectrum for specific and generic shapes of non-Gaussianities [34,35], to our knowledge, cosmological N-body simulations with non-Gaussian initial conditions specified by a bispectrum have been so far run only for local type of non-Gaussianity .
Given the importance N-body simulations have in modeling both the non-linear physics and observational effects that play such a crucial role in interpreting large-scale structure data, here we demonstrate how to set up and run N-body simulations with non-Gaussian initial conditions specified by a generic bispectrum.
This paper is organized as follows. In Sec. 2 we outline how to create a threedimensional non-Gaussian field with a given bispectrum. In Sec. 3 we explicitly work out the expressions to use in the four most popular non-Gaussian shapes discussed above and detail the steps for implementation. We test the non-Gaussian initial conditions in Sec. 4. In Sec. 5 we describe our simulations and use the local case as a benchmark for our implementation. In Sec. 6 we present our results. We conclude in Sec. 7.

Creating a 3D non-Gaussian field with a given bispectrum
The argument of [34] (see also [35]) 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 [38,39], here we follow a different route. Let's start from 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 Simulation were run e.g. for χ 2 initial conditions [36,37], but here we are concerned with inflationmotivated non-Gaussianities classified by their bispectrum shape.

If we define
with Φ * k the complex conjugate of Φ k , and use it in Eq. (2), we recover the identity of Eq. (1).
It is important to bear in mind (as already pointed out in [35]) that this procedure is not unique: in other words there may be more than one -non-equivalent-expression for Φ N G 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 [40,15,41]: 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: Note that Eq. (6) is not equivalent to Eq. (4)&(2). Both procedures yield Φ fields with the same bispectrum, but the Φ fields obtained are different. In fact it can be shown that Eq. (4)&(2) become equivalent to Eq. (6) only if B(k 1 , k 2 , k 3 ) −→ 6f NL P (k 2 )P (k 3 ) rather than using Eq. (7). This extra "degree of freedom" was used in Ref [35] 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. (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 [40,15,41] 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 [42,43,44,45]. These two shapes arise due to the interactions of the field driving inflation. On the other hand, modified initial state (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 orthogonal shape: 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: where 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: where F eq = − P (k 1 )P (k 2 ) + 2cyc. − 2[P (k 1 )P (k 2 )P (k 3 )] 2/3 (11) 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: where 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 if 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 Φ N G 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 [48] 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, Φ N G k , to the gravitational potential. For a generic, non-factorizable bispectrum, we can compute this field by disctretizing Eq. (4), i.e., for each grid point k in Fourier space, we loop over all grid points k : 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 Φ N G 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. 19, 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. (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. (5) and (6) where instead of the Φ field we have two auxiliary fields G and Q).
While for factorizable shapes the direct summation approach (Eq. 19) 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 [47] 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 ¶.

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. 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. (19) with the bispectrum given by Eq. (7) c) using again Eq. (19) 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 Fig. 4 and Fig. 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 ¶ Since we will be interested in the non-Gaussian to Gaussian ratio of our statistics, the detailed implementation of the initial displacement field (i.e., Zel'dovich or Second-Order Lagrangian perturbation theory) does not matter. 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%.
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 Φ N G k Φ N G −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 sec. V of [35]; 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- 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. 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 Fig. 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 Φ N G k Φ N G −k , which is further explored in the Appendix A. We observe 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 allat 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. [49]). In 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 Tab. 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 Sec. 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 [50]. Our choice of cosmology is a flat ΛCDM cosmology, which is consistent with the seven-year WMAP results [51]. 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. The results are presented in Fig. 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.7 and k 0.3 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 [52,53]. 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, we find that the ratio of the mass functions is not affected by the low mass resolution. This is demonstrated in Fig. 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, Fig. 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 Sec. 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., [27,28] and the bispectrum is left to forthcoming work [46].

Non-linear power spectrum
High precision in the theoretical prediction of the nonlinear matter power spectrum at k ∼ 1 h/Mpc and above is needed to derive unbiased results from the data of upcoming weak-lensing surveys e.g., [54]. 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 nonlinear power spectrum at the few percent level. In Fig. 6 (local), Fig. 7 (equilateral), and Fig. 8 (enfolded), we show the ratio of the non-Gaussian and Gaussian power spectrum for different type of non-Gaussianities with f NL = 500 (green), 250 (red), and 100 (blue) at redshift z = 1 (open squares) and z = 0 (filled circles). The black dashed lines correspond to the perturbation theory prediction using the timerenormalization group (TRG) approach [55]. For z = 1 and k 0.2 h/Mpc, the TRG    Figure 7. Same as in Fig. 6 but for the equilateral type of non-Gaussianity.  Figure 8. Same as in Fig. 6 but for the enfolded type of non-Gaussianity.
predictions agree very well with the results of the N-body simulations. At z = 0, perturbation theory slightly overpredicts the effect on scales around k ∼ 0.1 h/Mpc before it breaks down on smaller scales (k 0.2 h/Mpc).
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 [56].

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 [57]. Hence, galaxy cluster surveys offer in principle the possibility to probe primordial non-Gaussianities. In Fig. 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 N G (M ), provided in [20] and [18] into the corresponding ratio of the cumulative mass functions, R N G (> M ), we use the fitting formula of [58] for the SO halo mass function,  [20] and the red lines represent the model of [18].
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. 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 [59,60], 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. 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 lowmass end. However, for all shapes, the qualitative behavior is captured well by the analytic models. There is an indication that using spherical overdensities halos a value q ∼ 0.9 yield a slightly better fit to the 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 issue, 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.

Conclusions and discussion
We have addressed the issue of setting up generic non-Gaussian initial conditions for Nbody 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., [63,64,65] 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 [35].
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 halos need a q correction factor much closer to unity (if any correction at all) than Friends-of-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. Appendix 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 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 The term Φ G Φ N G 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.
Since the term Φ N G Φ N G is of second order and the potential is small, one could think that this term is negligible for reasonable values of f NL . However, our formula for Φ N G involves integrals which -depending on the bispectrum-are divergent. In order to explore this further, we derive the general expression for Φ N G Φ N G . After applying Eq. (4) for Φ N G , using the definition of the power spectrum, and integrating over the delta function we obtain As a first example let us consider the case in which the bispectrum of the local type (Eq. 8) is used for B in Eq. (4). In this case, we get three terms where the first and second term are very similar to the ones discussed in [61,62]; the second term is just a k-independent renormalization, which for reasonable values of f NL is negligible small because of the truncation of the integral due to the finite volume of the simulation box, while the first term gives rise to a k-dependent renormalization which results in a change in the slope of the power spectrum. However, for realistic values of f NL this change in the slope is very small [62]. The last term, proportional to P 2 (k), causes the largest modification of the power spectrum. Even for still allowed values of f NL this term gives significant contributions to the power spectrum on large scales (see Fig. A1).
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 equilateral case (see Eq. 11) 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 Φ N G 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. 17) this contributes substantially to the power spectrum and causes the deviations visible in Fig. 3.