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

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

[    [
###### Abstract

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 on 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 fudge factor (which was introduced before but still lacks convincing physical explanation) is calibrated on simulations. The halo bias remains therefore an extremely promising tool to probe primordial 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.

a]Christian Wagner a,b]and Licia Verde

Prepared for submission to JCAP

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

• ICCUB-IEEC, University of Barcelona, Barcelona 08028, Spain.

• ICREA, Institució Catalana de Recerca i Estudis Avançats.

Keywords: cosmological simulations, large scale structure, initial conditions

## 1 Introduction

Inflation, the standard theory for the origin of the primordial density fluctuations, predicts that the primordial density field can be described, to a good approximation, as a Gaussian random field. This nearly Gaussian nature of the primordial fluctuations is well established observationally (for current constraints see, e.g., [1]). However, small deviations from Gaussianity are predicted by inflation. The level and form of the 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, . 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 , gives the level of non-Gaussianity.

Since the Universe is statistically homogeneous and isotropic, the wave vectors , , and have to form a triangle. The shape of the triangle, for which the bispectrum peaks, is conventionally used to characterize the form of the non-Gaussianity. The local type of non-Gaussianity, , where is gravitational potential in the Gaussian case, yields a bispectrum that peaks for squeezed triangles (e.g., ). 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].

Other inflationary models which involve higher-derivative operators (see e.g. [13, 14, 15]) give rise to a bispectrum that has its maximum for equilateral triangles ().

Models in which the initial state of the inflaton is different from the 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 ().

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 assumed111For an approach without assuming a shape for the bispectrum, see [19]., constraints on the corresponding value can be derived from the observed fluctuations in the CMB. The current constraints for the most commonly used shapes are , , (95% CL) [1], and (68% CL) [20] . The Planck mission [21] is expected to measure with an accuracy of .

Primordial 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 , which is smaller than the scales accessible by CMB measurements. Intriguingly, several groups [30, 31, 32] derived estimates of from the observations of very massive galaxy clusters, which are as large as (see, however, [33]).

The clustering of halos on large scales offers another possibility to probe 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 , 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 , 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, (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 with an uncertainty of [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 non-linearities 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 simulations— the 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.

## 2 Non-Gaussian halo bias

The halo bias describes the relation between the halo density contrast, , and the underlying dark matter density contrast, . In Fourier space we can write the relation as , where the random variable models the stochasticity coming from shot noise and possibly from physical halo formation processes. For Gaussian initial conditions the bias is scale-independent on large scales, for [49, 50, 51]. Analytic predictions and results of N-body simulations show that the linear222The linear bias, , is called linear as it corresponds to the linear coefficient in the Taylor expansion of the halo density field: , where both density fields are smoothed on a sufficiently large scale. Gaussian bias, , increases with redshift and halo mass . In addition the halo bias may depend on environment and mass accretion history [52, 53].

Primordial non-Gaussianity of the local type induces a scale dependence of the halo bias on large scales, which scales as . This particular scaling was first found in 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 [45, 44, 46, 47]. Reasonable agreement between theory and simulations was found if small corrections, which we discuss in detail below, are taken into account. In a number of subsequent papers further theoretical modelling was developed by applying (univariate and multivariate) local biasing in the framework of perturbation theory [55, 56, 57, 58, 59, 27]. In summary the non-Gaussian bias can be modelled to leading order as

 Δb≡bNGM,z(k,fNL)−bGM,z(k)=Δb1;M,z(fNL)+fNL[bG1;M,z+Δb1;M,z(fNL)−1]qδcD(z)FM(k)MM(k), (2.1)

where is the difference in the linear bias, which is non-zero, since a non-vanishing changes the halo number density and hence the linear halo bias for a fixed halo mass . Mainly, this effect leads to a scale-independent shift relative to the Gaussian () linear bias [44, 59, 60]. Additionally, it makes the scale-dependent second term in above equation slightly non-linear in [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 , the linear Lagrangian bias , and the redshift-dependent collapse threshold , where is the spherical collapse threshold in an Einstein-de Sitter universe, is the growth function normalized to at high redshift, and can either be regarded as a fudge factor introduced to improve the fitting to 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 and . connects the density with the primordial potential. The linear matter density contrast smoothed on a mass scale is related to the primordial potential by the Poisson equation,

 δm;M,z(k)=23T(k)k2ΩmH20WM(k)D(z)Φ(k)≡MM(k)D(z)Φ(k). (2.2)

Here, and are the present day matter fraction and the Hubble constant, respectively. The transfer function models the physics until the end of recombination and is normalised to unity on large scales. The function is the Fourier transform of the top-hat filter with a smoothing scale given by the mass . As , we see that scales as on large scales.

The “form factor”, , depends on the shape of the non-Gaussianity through the bispectrum of the primordial potential , [35]:

 FM(k)=18π2σ2MfNL∫dk′k′2MM(k′)∫1−1dμMM(|k+k′|)BΦ(k,k′,|k+k′|)PΦ(k), (2.3)

where is defined by , is the primordial power spectrum, and denotes the variance of linear density fluctuations on the mass scale at redshiht . Considering the limit of , we see that on large scales the form factor depends primarily on the squeezed limit of the bispectrum. For the local type of non-Gaussianity, the “form factor” becomes scale and mass independent on large scales, . Hence, we obtain the aforementioned scale dependence of the halo bias in the case of local 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 , which was introduced in [44], is a relatively small contribution coming from the fractional difference in the Gaussian and non-Gaussian non-linear matter power spectrum. This effect becomes larger at smaller scales but remains below a few percent on all scales for observationally allowed values of .

Some additional effects are also missed in Eq. (2.1): The impact of the mass accretion history of the halos —in particular recent major mergers— on the 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, , changes the amplitude of the non-Gaussian correction by , where varies roughly between and from very young to very old halos [63]. For recent major mergers, takes the value [37, 63].

In addition, higher-order effects in the biasing description (e.g., terms of the order [55, 64, 59]) and in the primordial non-Gaussianity (e.g., non-Gaussianity at the cubic order quantified by [64, 59]) are neglected in Eq. (2.1).

Finally, general relativistic effects are not taken into account, which might become important on very large scales [65, 66, 67, 36, 68, 69].

So far we discussed the scale-dependent halo bias caused by the local type of non-Gaussianity with a scale-independent . The effect of a scale-dependent 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 . 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].

## 3 Physical shapes vs. templates

The functional form of the primordial bispectrum, , predicted by inflationary models can be quite complex. In particular, the bispectrum is usually not separable, i.e. it cannot be written as a finite sum of terms which are factorizable as a product of functions of , and . Separability is important for computationally efficient analysis and simulation of both CMB maps [13, 19] and 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 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 equilateral333This template peaks when the wavenumbers form a equilateral triangle. template, which was introduced in [74] as a separable approximation to the shape functions of DBI inflation [14], inflation with higher-derivative terms [13], and ghost inflation [15]:

 BeqlΦ(k1,k2,k3)=6feqlNL(−Floc−2FA+FB), (3.1)

with444The original expression of the equilateral template was proposed for a scale-invariant power spectrum and given in terms of the wavenumbers . Here we generalize it to non scale-invariant power spectra by substituting with with proper normalization. We use the same substitution ansatz in the case of other bispectra.

 Floc = PΦ(k1)PΦ(k2)+PΦ(k2)PΦ(k3)+PΦ(k1)PΦ(k3) (3.2) FA = [PΦ(k1)PΦ(k2)PΦ(k3)]2/3 (3.3) FB = P1/3Φ(k1)P2/3Φ(k2)PΦ(k3)+5perm. (3.4)

and the dependence of on the three wavenumbers is understood. The primordial power spectrum is specified by the amplitude at and the spectral index :

 PΦ(k)=A0(k/k0)ns−1k−3. (3.5)

The equilateral template was also adopted in other models. More generally, the form of the bispectrum is interesting because of its direct connection to the (self) interactions of the inflaton. In the effective field theory of inflation [75], there is a direct connection between the coefficients of operators of the Lagrangian for the inflationary fluctuations and the shapes of primordial non-Gaussianity. In [18], it was shown that the bispectrum of single-field inflation with an approximate shift symmetry protecting the Goldstone boson is a linear combination of two different shapes. The two shapes, and , which correspond to the interaction operators and (see [18] for details), are both close to the equilateral template. Especially, the cosine of with the equilateral template is very close to one. For this reason, we denote this shape also by “SSZ equilateral”.555Previously, in the framework of general single field inflation, [9] obtained two independent shapes which span the same space as the two shapes of [18].

We now turn to the question how well the equilateral template captures the behaviour of these physical shapes in the squeezed limit. Choosing and using the following notation, , , and , we obtain for the equilateral template

 BeqlΦsqueezed = 12feqlNL[P1/3Φ(k)P5/3Φ(k′)−(ns−43)2μ2PΦ(k)k2PΦ(k′)k′−2] (3.6) = 12A20feqlNL⎡⎢⎣⎛⎝k13k′53k20⎞⎠ns−1−(ns−43)2μ2(kkk20′)ns−1⎤⎥⎦k−1k′−5 (3.7) = 12A20feqlNL(1−μ2)k−1k′−5∼k−1, (3.8)

where in the last line a scale-invariant power spectrum was assumed, i.e. . We find a similar scaling of the bispectrum in the squeezed limit also for the physical shapes generated by the aforementioned inflationary models. However, the exact scaling and its amplitude differ slightly from case to case. In the case of a scale-invariant power spectrum, the scaling reduces for all models to with a model-dependent amplitude. For example for the bispectrum of DBI inflation (see [8] for an explicit expression) we get

 BDBIΦsqueezed=997A20feqlNL(1−511μ2)k−1k′−5, (3.9)

while the scaling of the bispectrum corresponding to [18] is

 BSSZ eql.Φsqueezed=21617A20feqlNL(1−58μ2)k−1k′−5. (3.10)

Here, we used the equilateral normalization convention for all shapes, i.e. .

In Fig. 1, the scale-dependent part of the non-Gaussian halo bias, , is shown for different bispectra of the equilateral type. The dashed lines depict 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 and , respectively, in agreement with the above expressions for the squeezed limit. The solid lines show computed with a power spectrum with . The upturn and the wiggles on smaller scales stem from the transfer function included in . Note, in contrast to the local type of non-Gaussianity, 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 and , respectively.

For other templates, as we discuss now, the situation becomes more problematic. The orthogonal666The orthogonal template was introduced as a separable approximation to the linear combination of the shapes and that is orthogonal (i.e. has a vanishing cosine) to template [18]

 BorthΦ(k1,k2,k3)=6forthNL(−3Floc−8FA+3FB), (3.11)

and the enfolded777This 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. . template [16]

 BenflΦ(k1,k2,k3)=6fenflNL(Floc+3FA−FB)=12(BeqlΦ−BorthΦ), (3.12)

scale both as in the squeezed limit, while the corresponding physical shapes scale as and , respectively. In Fig. 2 we show again the scale-dependent part of the non-Gaussian bias, , computed from these commonly used templates and the corresponding physical shapes. For comparison, the scale dependence obtained from the local type of non-Gaussianity, , 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 and , 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 is rescaled by . However, on scales of the different halo mass dependence lifts this degeneracy to some extent. The orthogonal shape of [18] (SSZ orth.) —having the same scaling in the squeezed limit as the equilateral shapes— does not lead to a 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 or as in the squeezed limit, while the scaling of the enfolded and orthogonal templates is .

In [19] a mode expansion of the bispectrum in separable functions was introduced. This expansion can be applied to arbitrary 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 in the squeezed limit, irrespectively of the scaling of the original bispectrum.

New separable templates which are not only overall good approximations to the original bispectrum but also faithfully capture the scaling in the squeezed limit can most probably be constructed (e.g [20]). However, we will see in the next section that separability of the bispectrum does not speed up our current implementation of generating the initial conditions for 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.

## 4 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

 ⟨Φk1Φk2Φk3⟩=(2π)3δD(k1+k2+k3)BΦ(k1,k2,k3). (4.1)

We shortly review our method. First, we decompose the potential field in Fourier space, , in a Gaussian and a non-Gaussian part

 Φk=ΦGk+ΦNGk. (4.2)

The Gaussian field is statistically completely described by the power spectrum

 ⟨ΦGk1ΦGk2⟩=(2π)3δD(k1+k2)PΦ(k1). (4.3)

As the non-Gaussian contribution is small, the leading order of the bispectrum is given by

 ⟨Φk1Φk2Φk3⟩=⟨ΦGk1ΦGk2ΦNGk3⟩+⟨ΦGk1ΦNGk2ΦGk3⟩+⟨ΦNGk1ΦGk2ΦGk3⟩. (4.4)

Using the following ansatz for

 ΦNGk = 16(2π)3∫d3k2d3k3BΦ(k,k2,k3)δD(k+k2+k3)Φ∗Gk2Φ∗Gk3PΦ(k2)PΦ(k3) (4.5) = 16(2π)3∫d3k2BΦ(k,k2,|k+k2|)Φ∗Gk2PΦ(k2)ΦGk+k2PΦ(|k+k2|),

with the complex conjugate of , one can easily show that

 ⟨ΦGk1ΦGk2ΦNGk3⟩=13(2π)3BΦ(k1,k2,k3)δD(k1+k2+k3). (4.6)

Finally, we recover at leading order Eq. (4.1) by virtue of the symmetry property of the bispectrum.

One problem of the above ansatz is that the non-Gaussian contributions to the power spectrum can become very large on large scales. 888On large scales, we require the non-Gaussian contributions to the power spectrum to be negligible (or renormalizable), in order to be in agreement with CMB measurements. Naively, one expects to be much smaller than the Gaussian contribution . However in the case of certain bispectra, diverges more strongly for small wavenumbers than the Gaussian power spectrum. With the above ansatz for the non-Gaussian contribution to the power spectrum reads

 ⟨ΦNGkΦNGq⟩=118δD(k+q)∫d3k′B2Φ(k,k′,|k+k′|)PΦ(k′)PΦ(|k+k′|). (4.7)

In practice, there is a minimum wavenumber , the largest mode which fits in the simulation box, and a maximum wavenumber , which corresponds to the smallest wavelength still resolved by the used grid size. These wavenumbers provide large and 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, , 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 , , and 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.999In [73] the mode expansion technique of [19] was used to develop a computationally efficient method to generate initial conditions for N-body simulations for a given bispectrum. However, due to the scaling in the squeezed limit of the expanded bispectrum, also in this case, non-Gaussian contributions to the power spectrum of the form arise on large scales.

For the local type of non-Gaussianity, we found in [48] that our ansatz becomes identical to the real space formulation , which does not lead to the aforementioned divergences, when we use in Eq. (4.5) instead of . Here we generalize this modification by writing it in the following way

 BΦ(k1,k2,k3)⟶BΦ(k1,k2,k3)3PΦ(k2)PΦ(k3)PΦ(k1)PΦ(k2)+PΦ(k2)PΦ(k3)+PΦ(k1)PΦ(k3). (4.8)

Inserting this expression in Eq. (4.5) yields the modified ansatz for (see also [72])

 ΦNGk = 12(2π)3∫d3k′BΦ(k,k′,|k+k′|)Φ∗Gk′ΦGk+k′PΦ(k)PΦ(k′)+PΦ(k′)PΦ(|k+k′|)+PΦ(k)PΦ(|k+k′|), (4.9)

which by using Eq. (4.4) reproduces Eq. (4.1).

Note that the integrand in the modified ansatz is in general not separable. This means that the integration cannot be written in a simple way as a sum of convolutions, which can be efficiently computed by Fast Fourier Transformations (FFTs). However, a possible avenue to still factorize our modified ansatz is the use of Schwinger parameters (see [76] for details). This is in contrast to our original ansatz, which is separable as long as the bispectrum is separable. However, the modified ansatz does not lead to contributions to the power spectrum, which diverge for , since they are suppressed by :101010Nevertheless, without cutoffs the integral can still be divergent, e.g. in the local case this would lead to a renormalization of the amplitude. In practice, the cutoffs given by the box size and grid size render this contribution negligible.

 ⟨ΦNGkΦNGq⟩=12δD(k+q)∫d3k′B2Φ(k,k′,|k+k′|)PΦ(k′)PΦ(|k+k′|)[PΦ(k)PΦ(k′)+PΦ(k′)PΦ(|k+k′|)+PΦ(k)PΦ(|k+k′|)]2. (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 into the linear density field with the help of the transfer function, which we compute with CAMB [77], and the Poisson equation:

 δk=23k2T(k)D(z)ΩmH0. (4.11)

Finally, the density is sampled with particles. The positions and velocities of the particles are derived from the displacement field at the initial redshift using the Zel’dovich Approximation or second-order Lagrangian perturbation theory [78, 79, 80].

## 5 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 and a particle load of one billion particles. This enables us to probe large scales () and resolve halos with masses above ( particles).

The downside of our modified ansatz for the initial non-Gaussian part of the potential, 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 mode the following sum has to be computed

 ΦNGk=12 ∑k′BΦ(k,k′,|k+k′|)Φ∗Gk′ΦGk+k′PΦ(k)PΦ(k′)+PΦ(k′)PΦ(|k+k′|)+PΦ(k)PΦ(|k+k′|), (5.1)

where runs over all grid points in Fourier space. If denotes the number of grid points in one dimension, then the computational cost of computing for all modes scales as . This scaling renders the usage of large grid sizes infeasible. Even for a moderately small grid size of , the generation of the 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, , which samples only the Gaussian contribution to the potential, and one smaller grid for the non-Gaussian part.111111To compute the non-Gaussian part, we use Eq. (5.1) where we only sum over modes which lie within the smaller grid. 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 . 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, .

We adopt a flat CDM cosmology with cosmological parameters consistent with current observational constraints [1], namely (the present-day matter fraction), (present-day baryon fraction), (the Hubble constant), (the spectral index of the primordial power spectrum), and (the rms of the linear density fluctuations at in spheres of ).

The simulations are carried out with Gadget-2 [81] taking only the gravitational interaction into account. The box size is , and particles with a softening length of are used to sample the density field. The simulations are started at an initial redshift of using 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 and of the halo and matter density field, respectively. This is achieved by first assigning the halos or particles to a regular grid () using the 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 and . We compute the halo bias by the ratio of the halo-matter cross power spectrum to the matter power spectrum

 b(k)=Phm(k)Pm(k), (5.2)

with and , where the average is taken over all modes which lie in a shell of radius and thickness . As expected, the clustering of halos with a mass below is affected by the smaller grid size used for the non-Gaussian . The grid spacing of the non-Gaussian field is . This corresponds approximately to the Lagrangian radius of halos of mass . The small difference in the halo clustering for more massive halos is consistent with shot noise. At higher redshift, the noise increases due the smaller number of halos per mass bin.

The other quantity of interest, which is also affected by the smaller 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 , 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, , above which the results are not affected by the low spatial resolution of the non-Gaussianity.

We conclude that, for halos more massive than , even the simulations which are set up with a small non-Gaussian grid () 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 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 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 . 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 .

## 6 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 -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.

### 6.1 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, , 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

 RNGMVJ(M,z)=exp(Δ3c(z)S3,M6σ2M)∣∣ ∣∣Δc(z)6√1−Δc(z)S3,M/3dS3,MdlnσM+√1−Δc(z)S3,M/3∣∣ ∣∣ (6.1)

and

 (6.2)

where is the linear mass variance on a mass scale at . The redshift dependence is modelled by the collapse threshold, , which includes a fudge factor, 121212This fudge factor was introduced in [46] and interpreted as a modification to the collapse threshold which also applies to the Gaussian mass function. Here, however, we do not adopt this interpretation and consider only as a fudge factor in the ratio of the non-Gaussian to the Gaussian mass function, which we calibrate using N-body simulations.. The non-Gaussianity is described by the normalized skewness of the linear density field, , on a mass scale at , . The normalized skewness scales linearly in . The magnitude and sign of depend on the shape of non-Gaussianity, while the mass dependence is only weakly dependent on the type of non-Gaussianity.

In Fig. 5, we show the fractional difference in the non-Gaussian and Gaussian mass functions derived from the simulations of the local type of non-Gaussianity. The four panels correspond to different redshifts, (top left), (top right), (bottom left), and (bottom right). The filled and open symbols represent the results for and , respectively. The corresponding analytic predictions of MVJ and LV (with the best-fit value for ) 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 , the data points do not follow the theory lines very closely. This is not of practical importance, as the difference is much smaller than observationally uncertainties in halo mass measurements. The applied fudge factor of for MVJ and for LV is in agreement with our findings in [48]. Other groups [46, 45], who worked with friends-of-friends (FOF) halos, used for both predictions, while [44], who applied a spherical-overdensity (SO) halo finder as we do, found for LV and for MVJ (V. Desjacques, private communication). For the standard definitions of FOF and SO halos, i.e. using a linking length of times the mean particle distance (FOF) and the 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 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 type are presented in Fig. 6. Here, the values are (filled symbols) and (open symbols). Note that these simulations used a rather small grid size () for the non-Gaussian initial density field. Hence, halos with masses below 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 to 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 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 . The third simulation, “eql. b)” (black triangles), used initial conditions generated with our original method, Eq. (4.5), with (see also App. A). Remember that in the case of the equilateral type, artificial contributions to the power spectrum are kept under control also when our original ansatz is used. As expected, the N-body results agree with each other for halo masses above . 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, 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 is in addition redshift dependent.

### 6.2 Halo bias

In order to explicitly explore the dependence on the halo mass, we divide the halos into five mass bins:

 3×1013M⊙/h< Mhalo <  6×1013M⊙/h 6×1013M⊙/h< Mhalo <1.2×1014M⊙/h 1.2×1014M⊙/h< Mhalo <2.4×1014M⊙/h 2.4×1014M⊙/h< Mhalo <4.8×1014M⊙/h 4.8×1014M⊙/h< Mhalo <9.6×1014M⊙/h.

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 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:

 Δb(k)=PNGhm(k)PNGm(k)−PGhm(k)PGm(k). (6.3)

As the non-Gaussian and the Gaussian simulation share the same realization of the initial Gaussian field, is almost free of sample variance and, in addition, has smaller shot noise than and individually. We estimate the error on directly from the distribution of in each bin (for details, see App. B).

Following Eq. 2.1, we model by

 Δb(k)=ΔbI+fNL[bG1+ΔbI−1]qδcD(z)FM(k)MM(k), (6.4)

where is the linear halo bias obtained from the Gaussian simulation on large scales (see App. B). The scale-independent shift, , 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]

 ΔbI=bNG1−bG1=−∂lnRNG(M)∂δc, (6.5)

with being the ratio of the non-Gaussian to the Gaussian mass function. Here, however, we treat as a free parameter and compare it later on with the prediction derived from the mass functions. We choose as the second free parameter. All other quantities in Eq. (6.4) are computed from the theory and are kept fixed.

In Fig. 8, we show as an example the effect of local non-Gaussianity on the halo bias for halos of mass at . Note that we plot . As is negative, this addition of is needed to still make use of the logarithmic scale.

The different line types visualize the effect of the different terms in Eq. (6.4). The solid lines show the best fit to the data (using all modes up to ) and include all terms given above. The short-dashed lines neglect appearing inside the square brackets in Eq. (6.4). The inclusion of this term makes the non-Gaussian bias non-linear in [59, 61], since depends on . The dot-dashed lines neglect completely. This scale-independent bias shift becomes important on smaller scales (), 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 , but were found in simulation outputs at . Consequently, the number of halos is smaller than in the previous figure and the residual shot noise is larger. The line types have the same meaning as before. On large scales, the halo bias scales as as predicted by the theory. Hence, with increasing wavenumber, the effect does not drop as rapidly as in the local case and extends to smaller scales.

Next, we present an example of the non-Gaussian halo bias induced by the equilateral shape. In Fig. 10, the halo bias corresponding to two different mass bins at is shown. As expected, the scale dependence is very weak and in agreement with the theoretical predictions. In particular, the observed mass dependence of the effect is consistent with the model predictions.

After having discussed for each type of 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 are presented. Different colours correspond to different redshifts: (red), (green), (blue), and (magenta). Triangles, boxes, and circles depict the three different realizations of the initial Gaussian random field used for the generation of the initial conditions. In the case of the local and orthogonal type, the open symbols correspond to and , respectively. Filled symbols show the results for (local), (orth.), and (eql.). For clarity, the points of each mass bin are spread over the range of each mass bin.

For the local type of non-Gaussianity, we recover within the error bars the same value, which was needed to bring the mass functions in agreement with analytic predictions. For the orthogonal shape, the effect is suppressed by a factor which is similar to the one we found for the non-Gaussian mass function, . 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 . 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 normalized by are shown. The colour and symbol coding is the same as before. As open and filled symbols (corresponding to different values) are consistent with each other, we can infer that the scale-independent shift is linear in for the values probed. The solid lines represent the predictions from Eq. (6.5) using the LV mass function ratio, , and taking the measured fudge factor from the mass function into account. Keeping in mind that the LV mass function is not in all cases a good fit to the mass functions derived from our simulations (see for example bottom right panel of Fig. 6), the agreement between the best-fit values and the predicted bias shift is very reasonable.

## 7 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 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 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 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 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 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 -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, , 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 shape-dependent normalization factor is included: effectively re-scales . 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 modifies the mass function and thus modifies also the linear bias. However, a Gaussian distribution with the same mass function would have a linear halo bias modified in the same way.

In most inflationary models considered so far, the bispectrum scales either as the equilateral shape () or as the local shape () in the squeezed limit (see, however, the 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 can be re-interpreted as limits on e.g., modified-initial-state type of non-Gaussianity, by an appropriate rescaling of : .

## Acknowledgements

We would like to thank Lotfi Boubekeur for useful discussions. The simulations and computations were performed on the cluster Hipatia (UB computing facilities and ERC grant FP7- IDEAS Phys.LSS 240117). We are grateful to the Centro de Ciencias de Benasque Pedro Pascual where the very last stages of the work were carried out. CW is supported by MICINN grant AYA2008-03531. LV acknowledges support from FP7-PEOPLE-2007-4-3-IRG n. 202182 and FP7- IDEAS Phys.LSS 240117 and MICINN grant AYA2008-03531.

## Appendix A Simulations of the equilateral type using two different approaches

In this section we compare the results of two N-body simulations with non-Gaussian initial conditions of the equilateral type. Both simulations have a box size of and a particle load of