# Transients from Initial Conditions in Cosmological Simulations

Martín Crocce, Sebastián Pueblas and Román Scoccimarro
Center for Cosmology and Particle Physics, Department of Physics, New York University, New York, NY 10003
###### Abstract

We study the impact of setting initial conditions in numerical simulations using the standard procedure based on the Zel’dovich approximation (ZA). As it is well known from perturbation theory, ZA initial conditions have incorrect second and higher-order growth and therefore excite long-lived transients in the evolution of the statistical properties of density and velocity fields. We also study the improvement brought by using more accurate initial conditions based on second-order Lagrangian perturbation theory (2LPT). We show that 2LPT initial conditions reduce transients significantly and thus are much more appropriate for numerical simulations devoted to precision cosmology. Using controlled numerical experiments with ZA and 2LPT initial conditions we show that simulations started at redshift using the ZA underestimate the power spectrum in the nonlinear regime by about at respectively, whereas the mass function of dark matter halos is underestimated by at () and at (). The clustering of halos is also affected to the few percent level at . These systematics effects are typically larger than statistical uncertainties in recent mass function and power spectrum fitting formulae extracted from numerical simulations. At large scales, the measured transients in higher-order correlations can be understood from first principle calculations based on perturbation theory.

###### keywords:
large-scale structure of the universe – methods: numerical.
pagerange: Transients from Initial Conditions in Cosmological SimulationsBpubyear: 2006

## 1 Introduction

The growth of structure in the universe is a fundamental aspect of cosmology, and can be used to constrain the nature of primordial fluctuations and the energy contents of the universe (e.g. dark matter and dark energy). Numerical simulations are a very important tool in making theoretical predictions of cosmological observables at small scales when fluctuations are nonlinear. Fluctuations are followed from their primordial seeds through the radiation era to the matter dominated era and current accelerating phase by using linear Boltzmann equation solvers (Peebles & Yu, 1970; Bond & Efstathiou, 1984; Ma & Bertschinger, 1995; Seljak & Zaldarriaga, 1996; Lewis et al., 2000), which are then coupled to N-body simulations to follow their nonlinear evolution. This coupling takes the form of final conditions predicted from linear Boltzmann codes input as initial conditions into N-body simulations at some redshift .

There are competing requirements on the best choice of . One would like to start a simulation as early as possible, where perturbations are closest to their linear value but this can be challenging from numerical reasons since one is trying to impose a small perturbation spectrum that can be overwhelmed by sources of noise such as numerical roundoff error or shot noise of the particles. Also, when running simulations that include a baryonic component, the initial conditions simplify considerably the lower is, in which case the baryons had more time to fall into the potential wells of the dark matter.

Despite the rapid progress in the development of different numerical N-body techniques for solving the nonlinear growth of structure in the universe (see e.g. Bertschinger (1998) for a review), the way initial conditions are input into these codes has essentially been the same since the first developments in cosmological simulations over twenty years ago (Klypin & Shandarin, 1983; Efstathiou et al., 1985), by using the Zel’dovich (1970) approximation (hereafter ZA).

As observations of the late universe quickly approach percent-level accuracy, it is vital to examine to what extent theoretical predictions can accurately go from primordial fluctuations to different stages of nonlinear evolution in recent epochs, e.g. redshifts . Studies of linear Boltzmann solvers (Seljak et al., 2003) have shown that independent codes agree to within 0.1%, whereas for N-body codes the agreement, while much more difficult to quantify, is of the order of few to ten percent depending on the statistics and scales considered (Splinter et al., 1998; Frenk et al., 1999; Heitmann et al., 2005).

The link between Boltzmann and N-body codes, or how to set initial conditions for nonlinear evolution, has been surprisingly little studied. The accuracy of the ZA has been essentially taken for granted and instead studies of initial conditions have concentrated on other (important) issues, related to how to best implement the ZA in different situations (Bertschinger, 1995; Salmon, 1996; Pen, 1997; Bertschinger, 2001; Yoshida et al., 2003; Sirko, 2005).

The rationale for this state of affairs seems to be that starting a simulation at “high enough” redshift with linear perturbation theory (which the ZA is, in Lagrangian space) should be safe. A rule of thumb for the choice of starting redshift often used is that the rms density fluctuations at the interparticle distance scale be of order . As a result, low resolution simulations are typically started much later than high resolution ones.

However, one of the lessons of perturbation theory (hereafter PT, see Bernardeau et al. (2001) for a review) is that the ZA has a very limited validity, reproducing only the correct linear growth and underestimating the skewness and higher-order moments (Grinstein & Wise, 1987; Juszkiewicz et al., 1993; Bernardeau, 1994; Catelan & Moscardini, 1994; Juszkiewicz et al., 1995; Baugh et al., 1995; Colombi et al., 1996), with deviations being particularly severe for velocity field perturbations because of the ZA failure to conserve momentum.

As a result of this, starting N-body simulations from ZA initial conditions leads to incorrect second and higher-order growing modes; this excites nonlinear decaying modes (called transients) that describe how long it takes for the correct dynamics to recover the actual statistical properties of density and velocity fields (Scoccimarro, 1998). This transient evolution is well understood at large scales using PT and describes the anomalous behavior seen in measurements of higher-order moments in N-body simulations started from ZA initial conditions (as shown in Scoccimarro (1998) and demonstrated in much more detail below).

The importance of transients is not just of interest for studies involving higher-order moments, as we show in this paper. Underestimation of non-Gaussianity means that the tails of the probability distribution of densities and velocities will be artificially suppressed, which implies that rare events will be even more rare in a simulation that has transients. Since typically large volumes (and thus typically low resolution) are used to study rare events, the dark matter halo mass function at the high-mass end will be particularly affected by transients given that low-resolution simulations are often started later than high-resolution ones. In addition, the transient behavior of large-scale higher-order moments implies that nonlinear couplings take a significant time to develop their full strength. These are the same couplings that act at small scales and give rise to the nonlinear evolution of the power spectrum, therefore transients should also affect the power spectrum at nonlinear scales.

In this paper we revisit the issue of transients from ZA initial conditions and study numerically the effects on higher-order moments, the mass function of dark matter halos and the power spectrum. As proposed in Scoccimarro (1998), transients can be reduced significantly by using second-order Lagrangian PT (hereafter 2LPT). This leads to correct second-order growth factors and greatly improves the behavior of higher-order moments. Furthermore, transients in this case decay as inverse square power of the scale factor, much faster than the inverse of the scale factor in the ZA case. We study the improvement brought by using 2LPT instead of ZA in the behavior of transients and quantify their magnitude in different statistics commonly used.

The paper is organized as follows. In section 2 we review how transients arise and present predictions that are contrasted with simulations later. Section 3 discusses the numerical simulations we run to quantify the transients in the nonlinear regime. In section 4 we present the main results of the paper on the magnitude of transients in higher-order cumulants, power spectrum, bispectrum, mass function and bias of dark matter halos. We summarize the results and conclude in section 5. In appendix A we discuss second-order Lagrangian perturbation theory and in appendix B we derive analytically the evolution of perturbations toward a glass configuration and discuss glass initial conditions from the point of view of transients.

## 2 Transients

In this section we discuss why transients exist and how they affect the statistics of density and velocity fields at large scales, following the treatment by Scoccimarro (1998), to which we refer the reader for more details. The analytic results summarized here also give a framework to understand the measurements in simulations we discuss below at small scales. Throughout this paper we assume Gaussian primordial fluctuations. For the most part we assume a flat universe with and , we mention the more general case at the end of section 2.1.

### 2.1 Equations of Motion

The equations of motion for the evolution of dark matter can be written in a compact way by introducing the two-component “vector”

 Ψa(k,η)≡(δ(k,η), −θ(k,η)/H), (1)

where the index selects the density or velocity components, with being the Fourier transform of the density contrast and similarly for the peculiar velocity divergence . is the conformal expansion rate with the cosmological scale factor and the conformal time. The time variable is defined from the scale factor by

 η≡lna(τ), (2)

corresponding to the number of e-folds of expansion. The equations of motion in Fourier space can then be written as (we henceforth use the convention that repeated Fourier arguments are integrated over)

 ∂ηΨa(k)+ΩabΨb(k)=γ(s)abc(k,k1,k2) Ψb(k1) Ψc(k2), (3)

where

 Ωab≡[0−1−3/21/2], (4)

describes the linear structure of the equations of motion. In particular, the first row describes the linearized continuity equation and the second row describes linearized momentum conservation, with the two entries being the Poisson equation () and the Hubble drag term, respectively. The symmetrized vertex matrix describes the non linear interactions between different Fourier modes and is given by

 γ(s)121(k,k1,k2) = δD(k−k1−k2) (k1+k2)⋅k12k21, (5) γ(s)222(k,k1,k2) = δD(k−k1−k2) |k1+k2|2(k1⋅k2)2k21k22,

, and is zero otherwise, denotes the Dirac delta distribution. The formal integral solution to Eq. (3) is given by (see Scoccimarro (2000) for a detailed derivation)

 Ψa(k,η) = gab(η) ϕb(k)+∫η0dη′ gab(η−η′) (7) ×γ(s)bcd(k,k1,k2)Ψc(k1,η′)Ψd(k2,η′),

where denotes the initial conditions, set when the linear growth factor and . The linear propagator is simply related to , and given by an inverse Laplace transform ()

 gab(η)=∮ds2πi esη (Ωab+sδab)−1 (8)

which given Eq. (4) yields

 gab(η)=eη5[3232]−e−3η/25[−223−3]. (9)

This describes the standard linear evolution of the density and velocity fields from their initial conditions specified by .

In the ZA, the Poisson equation is replaced by (Munshi & Starobinsky, 1994; Hui & Bertschinger, 1996), and thus Eq. (4) is replaced by,

 ΩZAab≡[0−10−1], (10)

which from Eq. (8) gives the linear propagator

 gZAab(η)=eη[0101]+[1−100]. (11)

The change in the linear propagator (a reflection of gravity being driven by the velocity divergence instead of density fluctuations) has important consequences, as its structure determines the generation of non-Gaussianity.

From Eq. (7) we can see why transients arise. For simplicity, let’s assume first that the equations of motion are linear, i.e. that and the second term in the right hand side of Eq. (7) is zero. Then, selecting growing mode initial conditions is very simple. All we have to do is to select , where is a Gaussian random field whose power spectrum follows from a Boltzmann solver and selects the growing mode upon contraction with in Eq. (9). Selecting the growing mode is easy since we can solve the linear evolution analytically (given by the first term in Eq. (7)), otherwise one would generally excite the decaying mode.

In the nonlinear case, we do not know an exact solution (that’s the purpose of running a simulation), therefore we cannot select initial conditions in the exact (nonlinear) growing mode. The best we can do is to try to use PT to select the fastest growing mode at each order in PT (2LPT achieves this to second-order, and does a good job for higher orders as well). Since the initial conditions will not be in the exact growing mode, decaying modes will be excited by nonlinear interactions in the second term of Eq. (7): these are what we call transients. This is analogous to what would happen in the linear case if was not the eigenvector , i.e. the decaying mode in the second term of Eq.(9) would be excited.

In the nonlinear case the problem of transients is worse than one might think based on linear theory because transients do not go away as fast as in linear theory, where the decaying mode is suppresed by compared to the growing mode. Indeed, if the ZA is used to set initial conditions transients are only suppressed by , so they take a long time to decay. Using better initial conditions, such as 2LPT, leads to suppression, a substantial improvement.

To describe transients in a quantitative way, we expand the initial and final conditions in terms of Gaussian linearly evolved density perturbations ,

 ϕa(k)=∞∑n=1ϕ(n)a(k),     Ψa(k,η)=∞∑n=1Ψ(n)a(k,η), (12)

with

 ϕ(n)a(k)=[δD]n I(n)a(k1,…,kn) δ0(k1)…δ0(kn), (13)
 Ψ(n)a(k,η)=[δD]n F(n)a(k1,…,kn;η) δ0(k1)…δ0(kn), (14)

where and . The kernels describe how the initial conditions are set (at ), e.g. if linear theory is used and for . The kernels describe how perturbations evolve after initial conditions are set, they obey the recursion relation (Scoccimarro, 1998),

 F(n)a(η) = n−1∑m=1∫η0dsgab(η−s)γ(s)bcdF(m)c(s)F(n−m)d(s) (15) + gab(η) I(n)b,

obtained by replacing Eq. (12) and Eqs. (13-14) into Eq. (7). For brevity we have suppressed the Fourier arguments. Note that as it should be, afterwards the fastest growing mode behaves as at each order . Transients determine how fast one goes from initial conditions given by to the regime where final conditions are in the fastest growing mode .

If the initial conditions are not set using linear theory, the kernels obey recursion relations that follow from the dynamics used to set them. For ZA initial conditions obeys,

 I(n)a = n−1∑m=1∫η−∞dsen(s−η)gZAab(η−s)γ(s)bcdI(m)cI(n−m)d,

with the condition . This can be understood from Eq. (15) by using linear theory initial conditions and keeping the fastest growing mode, formally done by evolving by an infinite amount of time (from ). Note that these kernels are independent of , that’s guaranteed due to the lower limit in the time integral. In the long time limit, the kernels approach , where obeys the recursion relation in Eq. (LABEL:recursionkernelsZA) with the linear propagator replaced by , Eq. (9). It is easy to check that if is replaced by such in Eq. (15), then is a solution at all times, i.e. there are no transients, as expected since one is setting initial conditions exactly, i.e. in the fastest growing mode.

Although most simulations are started by initial conditions decribed by the kernels in Eq. (LABEL:recursionkernelsZA) when using the ZA, some authors (Efstathiou et al., 1985; Power et al., 2003) define the velocities in a different way, by applying the ZA prescription to the perturbed density field after particles have been moved instead of the linear density field. In this case, it can be shown (Scoccimarro, 1998) that the velocity divergence kernels are not those in the ZA, but rather they are equal to the density kernels in the ZA. That leads to a reduction in the time scale of transients by a factor of about two (Scoccimarro, 1998). We will not consider this possibility here, as the standard ZA starts are most often used in the literature, and the improvement brought by using this prescription for the velocities does not fully address the limitations of the ZA (e.g. the second-order growth is incorrect, and the transients still evolve as the inverse of the scale factor).

As an example, the second-order density kernel reads from Eq. (15), with ,

 F(2)1(η) = e2η[57+x2(k1k2+k2k1)+27x2]+eη310(x2−1) (17) + e−3η/2335(x2−1),

which interpolates between the ZA value at and the exact value for large (that in between square brackets). Note that, as mentioned above, the dominant transient term is suppressed only by a single inverse power of the scale factor, so it takes a long time to disappear. If 2LPT is used to set initial conditions, the second-order kernels are reproduced exactly and there are no transients to this order. For higher-order kernels, 2LPT shows transients suppressed by inverse square powers of the scale factor (Scoccimarro, 1998), and as we see below they are greatly reduced compared to the ZA case.

In making predictions that evolve all the way to one must take into account the dependence on cosmological parameters neglected so far. A simple and rather accurate way of doing this is to redefine in Eq. (1), as the second component, where and is the growth factor with the exact cosmological dependence. When this is done, the only change in the equations of motion is in Eq. (4), and given that (Bouchet et al., 1995) for flat models with a cosmological constant this can be well approximated by the value of since during most of the time evolution is very close to unity. This leads to the same kernels as in the case but with the scale factor replaced everywhere by the growth factor .

### 2.2 Statistics

We now discuss quantitatively the impact of transients on statistics of the density field at large scales. Since linear growing modes are reproduced exactly, the power spectrum at large scales is not affected by transients. The effect of transients a large scales is only felt by higher-order spectra, such as the bispectrum, or cumulants, e.g. the skewness and kurtosis. However, as these quantities characterize non-linear couplings, they are useful for understanding the results at small scales where nonlinearities play a role in the power spectrum, which will then be affected by transients, as we shall see using simulations in section 4.2.

The power spectrum , bispectrum and trispectrum are given by,

 ⟨δ(k)δ(k′)⟩ = δD(k+k′) P(k), (18) ⟨δ(k1)δ(k2)δ(k3)⟩ = δD(k123) B123, (19) ⟨δ(k1)δ(k2)δ(k3)δ(k4)⟩c = δD(k1234) T1234, (20)

with , etc., and the subscript “c” denotes a connected contribution, whereas the skewness and kurtosis are given by,

 S3=⟨δ3⟩⟨δ2⟩2,          S4=⟨δ4⟩c⟨δ2⟩3. (21)

In second-order perturbation theory the bispectrum reads

 B123=2F(2)1(k1,k2;η) P0(k1)P0(k2)+ cyclic, (22)

where denotes the power spectrum of the linear fluctuations at the initial conditions (i.e. in linear theory), and cyclic permutations over are understood. The first three cumulants are given by,

 ⟨δ2⟩=∫P(k)W2(kR)d3k, (23)
 ⟨δ3⟩=∫B123W1W2W3δD(k123)d3k1d3k2d3k3, (24)
 ⟨δ4⟩c=∫T1234W1W2W3W4δD(k123)d3k1…d3k4 (25)

where represents a top hat filter in Fourier space at smoothing scale , and . From these expressions it follows the skewness and kurtosis including transients (Scoccimarro, 1998)

 S3 = {347−γ1}+1a[4−γ1]+1a(γ1−265)+1235a7/2,
 S4 = 607121323−62γ13+7γ313−2γ23−81635a+28γ15a+18475a2 (27) +1312245a7/2−8γ15a7/2−15044725a9/2+1921225a7,

where we have written the time dependence in terms of the scale factor and

 γp≡−dplnσ2(R)dlnpR, (28)

with the variance of density perturbations, Eq. (23).

These results have also been derived more recently using a different approach based on the steepest descent method for the density PDF (Valageas, 2002). In Eq. (LABEL:S3trans) we have explicitly written down the different contributions: the ZA contribution from initial conditions in square brackets, the exact value (reached only asymptotically as ) in curly brackets, and the remaining two transient terms (which cancel the exact value when ). If initial conditions were set using Eulerian linear theory, the term in square brackets would not be present.

As discussed above, to include the dependence of transients on cosmological parameters we will make the replacement in all comparisons of Eqs. (LABEL:S3trans-27) with numerical simulatons.

Figure 1 evaluates Eq. (LABEL:S3trans) in two cases, one where ZA is used to set the initial conditions (solid), and the other when Eulerian linear theory is used (and thus at ). We see that if linear theory were used, to reach the correct value of within some tolerance, it would take about four times more expansion away from the initial conditions compared to the ZA case. The difference in the two cases is purely due to the initial value of , i.e. the rate of convergence to the right answer is the same (as transients go away as in both cases). In contrast, using 2LPT initial conditions improves the initial values of the parameters dramatically over the ZA case (being exact for ), and also improves substantially the rate of decay of transients to .

Linear theory initial conditions are relevant for codes that do no use particles, e.g. in the case where baryons are described by fields in a grid instead of particles. In such cases one must keep in mind that the magnitude of the transients we demonstrate with simulations below are very likely to be much worse, as reflected by the difference between the curves in Fig. 1.

The dependence of transients on the shape of the power spectrum, in Eqs. (LABEL:S3trans-27) through , shows that transients from ZA initial conditions are somewhat reduced as the spectrum becomes steep, since the gravitational interactions are dominated by large scales and this leads to fewer deflections, i.e. more coherent displacement fields better approximated by straight-line trajectories given by the ZA. The dependence on spectrum is discussed in detail in Scoccimarro (1998). Here we concentrate on a fixed shape given by a concordance cosmological model as discussed in the next section.

The lesson from Fig. 1 is that nonlinear couplings (in this case, second order) take a significant time to develop their full strength if initial conditions are not set accurately. As expected from perturbation theory, these weaker nonlinear couplings give a transient effect in the nonlinear power spectrum. The effect of transients is also felt rather strongly on rare events such as the abundance of high-mass clusters, which probe the tail of the density probability distribution characterized by the parameters. We demonstrate these effects below using numerical simulations.

## 3 Simulations

In order to study the effects of transients at small scales we ran a set of simulations with different initial conditions (ZA and 2LPT), starting redshifts , box sizes , softening lengths and number of timesteps as described in Table LABEL:SimTable. All our simulations have particles with cosmological parameters , , , and , and were run using the Gadget2 code (Springel, 2005).

All the internal parameters to the Gadget2 code used to run the different numerical experiments are also given in Table LABEL:SimTable. In summary, we have two sets of initial conditions (ZA and 2LPT), each of them started at two different redshifts with the purpose of checking the magnitude of transients. The difference between ZA and 2LPT initial conditions is explained in Appendix A in terms of Lagrangian displacements. In our implementation, particles are started from a grid and all differentiations are done in Fourier space, i.e. through . Another possibility would be to start particles from a “glass”, in appendix B we briefly discuss how to avoid inducing extra initial skewness (and thus transients) due to glass initial conditions.

A simulation without transients should have statistical properties (measured at redshift ) independent of . Also, if generating initial conditions with the ZA was accurate enough, certainly 2LPT would be as well, since PT converges very well at high redshift. If this were the case, no appreciable difference would be seen at low redshift between ZA and 2LPT initial conditions.

However, for ZA initial conditions, we observe that statistics at low redshift are different than those from 2LPT initial conditions started at the same redshift (in this case ), and we also see a clear dependence on (comparing with simulations). Furthermore, the magnitude of transients observed at large scales agrees with the perturbative predictions presented in the previous section. In addition, the expected suppression of nonlinear couplings at large scales acts at small scales in a manner consistent with what one expects physically.

How do we know the 2LPT initial conditions simulations are correct? At large scales we checked the numerical simulations against the PT predictions, and we see almost no dependence of the low redshift results on the starting redshift . To make this latter point more obvious, we decided to present results for a very late start, (compared to our standard choice ). Note we are not advocating starting simulations at , this is only to prove that 2LPT initial conditions have very little transient effects. 2LPT initial conditions do have transient effects, although much smaller than ZA. Starting at redshift makes them large enough to be measurable at , which allows us to estimate the magnitude of transients in our reference runs (see section 5 for a discussion of this).

Most of the results presented here are obtained by using the “high-resolution” runs (labeled HR in Table LABEL:SimTable), corresponding to . The “low-resolution” (labeled LR) simulations with are used to study the high-mass tail of the dark matter halo mass function in Section 4.5. In addition, we have studied if the magnitude of transients depends on the softening length by reducing the softening length by a factor of three in the LR runs (denoted by LRs), and we have also run HR simulations with twice smaller softening length and stricter time-integration and force accuracy parameters, labeled “high-quality” (HQ) runs in Table LABEL:SimTable. From these sanity checks we conclude that the magnitude of transients observed from ZA initial conditions we present below is insensitive to the choices we have made in the less time-consuming runs.

The runs with 2LPT initial conditions with are our most accurate simulations from the point of view of transients. We denote such realizations as reference runs. In order to quantify the magnitude of transients, we measure statistics at low redshift and compare them to the same measurements in these reference runs.

## 4 Results

### 4.1 Higher-Order Cumulants

We first start by comparing simulations for measurements of higher-order cumulants at large scales, skewness and kurtosis , which are sensitive to transients and can be understood analytically using perturbation theory, as presented in Section 2.2.

Figure 2 shows the measurement of and as a function of smoothing scale for the reference runs HR2LPT49, after averaging over 20 realizations. The agreement at large scales with the perturbative prediction (solid lines) is excellent, better than 1% for for . Note that all measurements at different redshifts lie on top of each other in the large-scale limit as predicted by tree-level PT (Juszkiewicz et al., 1993; Bernardeau, 1994), and depart from it at increasingly large scales as decreases as expected from one-loop PT (Scoccimarro & Frieman, 1996; Fosalba & Gaztañaga, 1998). Previous measurements in large-volume simulations such as the Hubble volume (Colombi et al., 2000) could not achieve such a clean measure of higher-order cumulants independent of due to transients effects from ZA initial conditions, as will become clear shortly.

Figure 3 shows what happens with 2LPT initial conditions when the starting redshift is changed to a very late start, . As expected theoretically, at large scales there should be no change whatsoever in since 2LPT reproduces exactly the second-order growing modes, and that’s indeed what is seen in the bottom panel. At third-order, 2LPT does underestimate slightly, and that manifests in a small transient behavior at large scales (top panel). Note that since third- and higher-order couplings are slightly smaller in 2LPT, one expects the loop corrections to and to be smaller than in the reference runs, that’s indeed what is seen in both panels in Fig. 3 at small scales: there is a slight underestimate of the small-scale cumulants, although by the time one reaches , that is well under 1%. Overall, the behavior of 2LPT initial conditions is remarkably stable given the very late start, . In fact, we picked such a late start to make the effects visible at all at .

Figures 4 and 5 show the corresponding situation with ZA initial conditions for , respectively. It is obvious from these figures that the ZA initial conditions runs are not stable at large scales, in particular, and ZA runs do not give the same results, with significant deviations at the few percent level. These deviations are in excellent agreement with the PT predictions, Eq. (LABEL:S3trans) for the skewness (bottom panels) and Eq. (27) for the kurtosis (top panels). From this we conclude that the simulations started with ZA initial conditions have transients which are well understood at large scales from first principles. Figs. 4-5 check the PT transients predictions more accurately than done before (Scoccimarro, 1998). We now discuss how these effects impact on clustering at small scales.

### 4.2 Power Spectrum

The power spectrum is the most widely used statistic. Nonlinear corrections to the power spectrum are controlled by the same nonlinear couplings that determine higher-order statistics at large scales, and thus we expect to see transient behavior in the nonlinear power spectrum as well.

Figure 6 shows that this is indeed the case. First, in each panel the top line denotes the 2LPT initial conditions runs, showing that in the 2LPT case there is almost no transient behavior, as expected from the results on the parameters at large scales. On the other hand, the ZA initial conditions runs show significant transient effects, up to at , which can potentially lead to serious systematic errors for precision studies of the high-redshift universe. For example, a ZA start at was used for the LCDM simulation from which the latest fitting formula for the nonlinear power spectrum (Smith et al., 2003) was derived. Our reference runs show a power spectrum for that is as much as higher at and as much as higher at when compared to the Smith et al. (2003) fitting formula. These deviations are in reasonable agreement with the expectations based on Fig. 6. A detailed comparison of the nonlinear power spectrum in our simulations to fitting formulae, one-loop PT and renormalized PT (Crocce & Scoccimarro, 2006) will be discussed elsewhere.

It is also important to note that a relatively high-redshift ZA start such as , which we ran as a test, leads only to mild improvements compared to ZA in Fig. 6, i.e. the suppression in the nonlinear regime is still of order , consistent with the slow scaling of transients in the ZA. We have checked that these results do not depend on the accuracy of the force, time integration, and softening. For example, using the HQ runs (see table LABEL:SimTable) for we obtain very similar results to those in Fig. 6, specifically within about for and for . Also, as shown in Fig. 6, the magnitude of the measured transients in the ZA and ZA runs are in reasonable agreement with their estimation by one-loop PT using the kernels in Eq. (15), particularly at high redshift where the agreement should be best.

Figure 7 presents results in redshift space for the ZA initial conditions runs. The combined effect of density and velocity transients leads to nontrivial behavior. As shown by Scoccimarro (1998) transients in the velocity field are larger than in the density field. For modes parallel to the line of sight (triangles in each of the panels in Fig. 7), small-scale velocities suppress power due to velocity dispersion. Therefore, although transients in real space make the power smaller, when mapped into redshift-space with velocities that are even more suppressed by transients than the density, the overall effect is to increase the power spectrum along the line of sight compared to the reference runs that have no transients. Since there is cancellation between density and velocity transients, the overall effect is somewhat weaker than in the real space power.

For modes perpendicular to the line of sight, on the other hand, the effect is opposite since these modes are unaffected by redshift-space distortions and show the same suppression in power seen in Fig. 6. This means that when the redshift space power spectrum is averaged over angles with respect to the line of sight to obtain the monopole, these opposite behaviors lead to cancellation that makes the monopole less sensitive to transients (see circles in Fig. 7). However, the anisotropy of the redshift-space power spectrum is affected by transients much more than the monopole.

### 4.3 Bispectrum

The bispectrum is also affected by transients, and from Eqs. (17) and (22) we see that at the largest scales where second-order Eulerian PT holds it should be affected most for triangles that are most different from elongated triangles (which correspond to and show no transient behavior). Figure 8 shows the bispectrum at for different initial conditions divided by that in the reference runs for triangles with and as a function of the angle between and . From this we see again that the 2LPT initial conditions show very little () transient effects, whereas the ZA initial conditions runs show the expected transient signature, minimal at elongated triangles and maximal at isosceles triangles. These can lead to systematic errors in the determination of bias and cosmological parameters, given present expected observational errors (Sefusatti et al., 2005).

### 4.4 The PDF

The suppression of higher-order moments by transients means that the probability distribution function (PDF) of density fluctuations is also affected. Figure 9 shows a comparison of the PDF of our reference run (2LPT ) and ZA initial conditions at redshift for a smoothing scale . Note how the non-Gaussian features of the PDF are changed by transients: the right tail is suppressed, whereas the left cutoff is enhanced, leading to smaller skewness. Figure 10 shows the ratio of the PDF for different initial conditions to that of the reference run, to better appreciate the differences. Although the differences at the left tail are more significant, the PDF is falling very steeply, thus in practice the differences seen at the right tail are more important. These correspond to high-density regions, and thus we expect these differences to impact the high-mass tail of the mass function of dark matter halos.

### 4.5 The Mass Function

Figures 11 and 12 show the ratio of mass functions for different initial conditions to the reference runs, where the LR simulations have been used at the high-mass end to probe halos with masses larger than . All masses have been corrected by the procedure described in Warren et al. (2005) to take into account the finite number of particles making up halos. We use the friends-of-friends algorithm with linking length parameter equal to to identify halos.

The effect in the suppression of the high-mass tail is clearly seen in the simulations that were started with ZA initial conditions. The mass function in our reference runs is 15% larger than the Warren et al. (2005) fitting formula for at (and thus about larger than the Sheth & Tormen (1999) mass function), whereas for is consistent with the Warren et al. (2005) fit. The Jenkins et al. (2001) fit is very similar to the Warren et al. (2005) fit in the mass range covered by our simulations. These differences, in good agreement with the deviations seen in Fig. 11, are thus most likely due to transients in the simulations used to derived these fitting formulae at the high-mass end, which were started from ZA initial conditions.

Such a suppression of the mass function is clearly important for current and future observations that use the abundance of massive clusters to probe cosmology. It is also worth noting that the suppression of the mass function high-mass tail due to transients can be made worse by requiring that simulations start at a redshift such that density fluctuations at the interparticle distance be a fixed number (often chosen to be about ). This means that the large volume simulations that are required to probe rare events will be started later and thus will be subject to more suppression due to transients.

Using the LRs simulations with smaller softening length (see table LABEL:SimTable) for we obtain very similar results to those in Fig. 11 for , specifically the average deviation in the eight bins with (excluding the last one that has substantial errors) is within of that in Fig. 11. Therefore the suppression seen in the mass function is robust to a rather large change in softening length.

Reed et al. (2003) studied the effects of the starting redshift on the mass function measured at high redshifts from numerical simulations initialized with the ZA. Their statistical errors are rather large, but they do seem to see the effects we report here when comparing simulations with . They conclude is safe for measurements at and find a suppression of the high-mass tail from the Sheth & Tormen (1999) mass function which becomes stronger with . Extrapolating our results to their range, we would expect an important contribution from transients for those choices of going exactly in the same direction with increasing . Therefore we suggest that high-redshift studies of the high-mass tail of the mass function should be carefully checked for transients.

### 4.6 Halo Bias

Finally, Fig. 13 shows how the clustering of dark matter halos is affected by transients at . The panels show the halo-halo correlation function for different initial conditions compared to those in the reference runs, for three different mass bins. Here we have used HR simulations for the first bin in mass, and the LR simulations for the two largest mass bins. Note that as halo mass becomes larger, corresponding to rarer events, the halo bias in the simulations with more transients is larger. This is consistent with the fact that transients suppress the abundance of rare halos making them more rare and thus more biased, as expected from the peak-background split.

## 5 Summary and Conclusions

In this paper we considered the use of different approximations to set up initial conditions for cosmological simulations, concentrating on typical initial redshifts , and studied their influence on low-redshift () statistics. The standard procedure in the literature is to use the Zel’dovich (1970) approximation (ZA), a linear solution in Lagrangian space, or even linear perturbation theory in Eulerian space (typically when simulating baryons using methods which are not based on particles).

At high redshifts () the naive expectation is that linear perturbation theories should be accurate enough, but they are not. The reason for this is that at high redshifts nonlinear corrections are relatively much more important than at low redshift because the relevant scales involved have a very steep spectrum (with the effective spectral index very close to ). An example of this is shown in Fig. 9, where the density PDF smoothed on at shows strong non-Gaussian features, even though naively these scales are in the linear regime (since the linear rms density fluctuation is only ). Modeling these nonlinear corrections using the ZA leads to a significant underestimate of nonlinear couplings and thus non-Gaussianity. This leads to a suppression of the density PDF tail for large densities that translates into a suppression of the high-mass tail of the mass function. The weaker nonlinear couplings lead to a smaller power spectrum in the nonlinear regime.

The mismatch between the approximation used to set initial conditions and the correct dynamics leads to transients, excitation of decaying modes that are long-lived and survive until low redshift, systematically biasing low the nonlinear couplings even after evolving with the correct dynamics of an N-body code for as much as 50-100 scale factors. Transients cannot be avoided, but they can be significantly reduced by imposing initial conditions with better approximations to the nonlinear equations of motion.

More specifically, in this paper we have shown that,

1. The measurements of statistics at low redshift for simulations started with ZA initial conditions depend significantly on the starting redshift , signaling the presence of transients.

2. This can be avoided to a large extent by using 2LPT initial conditions, which as we have demonstrated show very little dependence on even for extremely late starts such as . Given that such runs show at most 1% transient effects and that transients in 2LPT decay as inverse square of the scale factor, we can estimate that our reference runs with are accurate (in terms of transients) to better than a tenth of a percent for .

3. The simplest way to quantify transients is to follow how non-Gaussianity develops at large scales as nonlinear couplings reach their full strength during the transient evolution. Our results for higher-order cumulants at large scales confirm previous PT results (Scoccimarro, 1998; Valageas, 2002) with much better accuracy than done before.

4. The small scale behavior of transients is consistent with the reduction of nonlinear couplings observed at large scales, and we have checked it is not affected by the choice of softening length, force and time integration accuracy.

5. The use of the ZA at leads to suppressions in the power spectrum of order a few () to almost ten percent () in the nonlinear regime. These results are in reasonable agreement with an estimation of transients in the power spectrum based on one-loop PT.

6. ZA initial conditions also lead to a suppression in the high-mass tail of the mass function of dark matter halos. For this suppression is about 10% for at and at . The clustering of halos is also affected at the few-percent level.

An obvious alternative to using 2LPT initial conditions is to start with the ZA much earlier than . However, in order to compete with the improvement brought by using 2LPT very high redshifts are required. For example, Fig. 6 shows that a 2LPT start has approximately a 0.5% transient effect on the nonlinear power spectrum at , whereas a ZA start has a 2.5% effect. In order for the ZA to reduce the transients to the level of the 2LPT start, a ZA start would be required. Reaching the level of a 2LPT start, corresponding to our reference runs, will require about another factor of 16 in redshift (given that transients scale as in 2LPT compared to in the ZA). Other issues make starting simulations at very high redshifts a bit more challenging, e.g. depending on the N-body code algorithm, numerical errors may be more difficult to avoid at very high redshifts where perturbations are so small. In addition, initial conditions for baryons at such high redshifts must be carefully included as they are still falling into dark matter potential wells.

Although we have not studied this in any detail, simulations started with linear theory instead of the ZA are expected to be much more significantly affected by transients, e.g. a simulation started with linear theory at should be roughly equivalent to one started with the ZA at (see Fig. 1). Note that linear theory initial conditions are often used in codes that rely on grid methods with gas dynamics solvers. We expect such methods to show differences with e.g. SPH methods that use ZA initial conditons. When the dynamical variables are densities and velocities in a mesh (instead of information carried by particles), a way to reduce transients would be to use higher-order Eulerian perturbation theory.

While in this paper we have concentrated on statistics of the dark matter density perturbations, it would be interesting to know how much do transients affect other high-redshift probes, where transients are most likely to play an important role. In this regard, the most pressing aspect would be to check the impact of transients in the Lyman- forest power spectrum, because of the high accuracy demanded from simulations to obtain cosmological constraints from current data sets (Viel et al., 2006; Seljak et al., 2006).

The use of better initial conditions should be a useful tool in reaching future goals of simulating the nonlinear power spectrum to an accuracy of better than 1% percent needed for the next-generation of weak gravitational lensing surveys (Huterer & Takada, 2005).

The code we used to generate ZA & 2LPT initial conditions in this paper is publically available at http://cosmo.nyu.edu/roman/2LPT.

## Acknowledgments

We thank S. Colombi, H. Couchman, S. Habib, K. Heitmann, L. Hui, A. Jenkins, A. Klypin, A. Kravtsov, U. Seljak, R. Smith, V. Springel, R. Thacker, M. Warren, D. Weinberg, for useful discussions and E. Sefusatti, M. Masjedi, and R. Smith for comments on the manuscript. We thank NYU Information Technology Services for making its High Performance Computation Cluster available to us, and thank Joseph Hargitai for his technical assistance. We acknowledge use of the University of Washington N-body shop friends-of-friends code to find dark matter halos.

## References

• Baugh et al. (1995) Baugh, C.M., Gaztañaga, E., Efstathiou, G., 1995, MNRAS, 274, 1049
• Bernardeau (1994) Bernardeau, F., 1994, ApJ, 433, 1
• Bernardeau et al. (2001) Bernardeau, F., Colombi, S., Gaztañaga, E., Scoccimarro, R., 2001, Phys. Rep., 367, 1
• Bond & Efstathiou (1984) Bond, J.R., Efstathiou, G., 1984, ApJ, 285, L45
• Bertschinger (1995) Bertschinger, E., 1995, preprint (arXiv:astro-ph/9506070)
• Bertschinger (1998) Bertschinger, E., 1998, ARAA, 36, 599
• Bertschinger (2001) Bertschinger, E., 2001, ApJS, 137, 1
• Bouchet et al. (1995) Bouchet, F.R., Colombi, S., Hivon, E., Juszkiewicz, R., 1995, A&A, 296, 575
• Buchert (1994) Buchert, T., 1994, MNRAS, 267, 811
• Buchert et al. (1994) Buchert, T., Melott, A.L., Weiss, A.G., 1994, A&A, 288, 349
• Carlberg & Couchman (1989) Carlberg, R.G., Couchman, H.M.P., 1989, ApJ, 340, 47
• Catelan & Moscardini (1994) Catelan, P., Moscardini, L., 1994, ApJ, 426, 14
• Catelan (1995) Catelan, P., 1995, MNRAS, 276, 115
• Colombi et al. (1996) Colombi, S., Bouchet, F.R., Hernquist, L., 1996, ApJ, 465, 14
• Colombi et al. (2000) Colombi, S., Szapudi, I., Jenkins, A., Colberg, J., 2000, MNRAS, 313, 711
• Couchman et al. (1995) Couchman, H.M.P., Thomas, P.A., Pearce, F.R., 1995, ApJ, 452, 797
• Crocce & Scoccimarro (2006) Crocce, M., Scocimarro, R., 2006, PRD, 063519
• Efstathiou et al. (1985) Efstathiou, G., Davis, M., Frenk, C.S., White, S.D.M., 1985, ApJS, 57, 241
• Fosalba & Gaztañaga (1998) Fosalba, P., Gaztañaga, E., 1998, MNRAS, 301, 503
• Frenk et al. (1999) Frenk, C.S., White, S.D.M., Bode, P., Bond, J.R., Bryan, G.L., Cen, R., Couchman, H.M.P., Evrard, A.E., Gnedin, N., Jenkins, A., Khokhlov, A.M., Klypin, A., Navarro, J.F., Norman, M.L., Ostriker, J.P., Owen, J.M., Pearce, F.R., Pen, U.-L., Steinmetz, M., Thomas, P.A., Villumsen, J.V., Wadsley, J.W., Warren, M.S., Xu, G., Yepes, G., 1999, ApJ, 525, 554
• Grinstein & Wise (1987) Grinstein, B., Wise, M., 1987, ApJ, 320, 448
• Heitmann et al. (2005) Heitmann, K., Ricker, P.M., Warren, M.S., Habib, S., 2005, ApJS, 160, 28
• Hui & Bertschinger (1996) Hui, L., Bertschinger, E., 1996, ApJ, 471, 1
• Huterer & Takada (2005) Huterer, D., Takada, M., 2005, Astropart. Phys., 23, 369
• Jenkins et al. (2001) Jenkins, A., Frenk, C.S., White, S.D.M., Colberg, J.M., Cole, S., Evrard, A.E., Couchman, H.M.P., Yoshida, N., 2001, MNRAS, 321, 372
• Juszkiewicz et al. (1993) Juszkiewicz, R., Bouchet, F.R., Colombi, S., 1993, ApJ, 412, L9
• Juszkiewicz et al. (1995) Juszkiewicz, R., Weinberg, D.H., Amsterdamski, P., Chodorowski, M., Bouchet, F.R., 1993, ApJ, 442, 39
• Klypin & Shandarin (1983) Klypin, A.A., Shandarin, S.F., 1983, MNRAS, 204, 891
• Lewis et al. (2000) Lewis, A., Challinor, A., Lasenby, A., 2000, ApJ, 538, 473
• Ma & Bertschinger (1995) Ma, C.-P., Bertschinger, E., 1995, ApJ, 455, 7
• Munshi & Starobinsky (1994) Munshi, D., Starobinsky, A.A., 1994, ApJ, 428, 433
• Pen (1997) Pen, U.-L., 1997, ApJ, 490, L127
• Peebles & Yu (1970) Peebles, P.J.E., Yu, J.T., 1970, ApJ, 162, 815
• Power et al. (2003) Power, C., Navarro, J.F., Jenkins, A., Frenk, C.S., White, S.D.M., Springel, V., Stadel, J., Quinn, T., 2003, MNRAS, 338, 14
• Reed et al. (2003) Reed, D., Gardner, J., Quinn, T., Stadel, J., Fardal, M., Lake, G., Governato, F., 2003, MNRAS, 346, 565
• Salmon (1996) Salmon, J., 1996, ApJ, 460, 59
• Scoccimarro & Frieman (1996) Scoccimarro, R., Frieman, J.A., 1996, ApJS, 105, 37
• Scoccimarro (1998) Scoccimarro, R., 1998, MNRAS, 299, 1097
• Scoccimarro (2000) Scoccimarro, R., 2000, Ann. N.Y. Acad. Sci. 927, 13
• Scoccimarro (2000b) Scoccimarro, R., 2000b, ApJ, 544, 597
• Sefusatti et al. (2005) Sefusatti E., Crocce, M., Pueblas, S., Scoccimarro, R., 2006, preprint (arXiv:astro-ph/0604505)
• Seljak & Zaldarriaga (1996) Seljak, U., Zaldarriaga, M., 1996, ApJ, 469, 437
• Seljak et al. (2003) Seljak, U., Sugiyama , N., White, M., Zaldarriaga, M., 2003, PRD, 68, 083507
• Seljak et al. (2006) Seljak, U., Slozar, A., McDonald, P., 2006 preprint (arXiv:astro-ph/0604335)
• Sheth & Tormen (1999) Sheth R.K., Tormen, G., 1999, MNRAS, 308, 119
• Sirko (2005) Sirko, E., 2005, ApJ, 634, 728
• Smith et al. (2003) Smith R.E., Peacock, J.A., Jenkins, A., White, S.D.M., Frenk, C.S., Pearce, F.R., Thomas, P.A., Efstathiou, G., Couchman, H.M.P., 2003, MNRAS, 341, 1311
• Splinter et al. (1998) Splinter R.J., Melott, A.L., Shandarin, S.F., Suto, Y., 1998, ApJ, 497, 38
• Springel (2005) Springel V., 2005, MNRAS, 364, 1105
• Yoshida et al. (2003) Yoshida N., Sugiyama, N., Hernquist, L., 2003, MNRAS, 344, 481
• Warren et al. (2005) Warrren M.S., Abazajian, K., Holz, D.E., Teodoro, L., 2005, preprint (arXiv:astro-ph/0506395)
• White (1996) White, S.D.M., 1996, in Schaefer R., Silk J., Spiro M., Zinn-Justin J., eds., Cosmology and Large-Scale Structure. Elsevier, Amsterdam (arXiv:astro-ph/9410043)
• Valageas (2002) Valageas P., 2002, A&A, 385, 761
• Viel et al. (2006) Viel, M., Haehnelt, M.G., Lewis, A., 2006 preprint (arXiv:astro-ph/0604310)
• Zel’dovich (1970) Zel’dovich Ya.B., 1970, A&A, 5, 84

## Appendix A Second-Order Lagrangian Perturbation Theory (2LPT)

Here we present the basics of 2LPT, we refer the reader to the literature (Buchert, 1994; Buchert et al., 1994; Bouchet et al., 1995; Catelan, 1995; Scoccimarro, 2000b) for more information. A detailed step by step implementation of 2LPT initial conditions for numerical simulations is presented in Appendix D2 of Scoccimarro (1998). Here we assume and . For precision work one must take into account the cosmological dependence, including radiation as well.

In Lagrangian dynamics, particle positions are described by a displacement field so that,

 x=q+Ψ, (29)

where particle trajectories obey the equation of motion,

 d2xdτ2+H(τ)dxdτ=−∇Φ, (30)

where is the gravitational potential. This leads to the following equation of motion for the displacement field,

 J(q,τ)∇⋅[d2xdτ2+H(τ)dxdτ]=32H2(J−1), (31)

where we have used the Poisson equation and where is the Jacobian of the mapping,

 J(q,τ)=Det(δij+Ψi,j), (32)

and .

Solving Eq. (31) perturbatively leads to positions and velocities,

 x = q−a∇qϕ(1)−37a2∇qϕ(2) (33) v = −aH∇qϕ(1)−67a2H∇qϕ(2), (34)

where the potentials satisfy the Poisson equations (Buchert et al., 1994),

 ∇2qϕ(1)(q) = δ(q) (35) ∇2qϕ(2)(q) = ∑i>j{ϕ(1)ii(q)ϕ(1)jj(q)−[ϕ(1)ij(q)]2}. (36)

Setting in Equations (33-34) leads to the Zel’dovich (1970) approximation.

Going beyond 2LPT to third-order (3LPT) becomes more costly due to the need to solve three additional Poisson equations (Buchert et al., 1994; Catelan, 1995), and the improvement in higher-order statistics is modest (Scoccimarro, 2000b), leading to an improvement by only a factor of two in the scale factor needed to suppress transients (Scoccimarro, 1998), whereas going from ZA to 2LPT the improvement is more than one order of magnitude. However, 3LPT does provide a better behavior in underdense regions (Bouchet et al., 1995), and may be worth considering for studies sensitive to the statistical properties of voids.

## Appendix B On Glass initial conditions

In this paper, our numerical simulations were started from a grid. An alternative to a grid start is to displace particles from a configuration known in the literature as a glass, which is obtained starting from a Poisson distribution and running an N-body code with the sign of the acceleration reversed (Carlberg & Couchman, 1989; Baugh et al., 1995; Couchman et al., 1995; White, 1996; Smith et al., 2003).

In this appendix we derive analytically how perturbations evolve towards a glass, and what steps must be taken to avoid exacerbating the problem of transients. Reversing accelerations means that the equations of motion are still given by Eq. (3) but with

 Ωab≡[0−13/21/2], (37)

instead of Eq. (4). As a result the linear propagator, Eq. (8), reads

 gab(η) = e−η/4 {cos(√234η)[1001] (38) +1√23sin(√234η)[14−6−1]},

instead of Eq. (9). The “growth” factor for density perturbations can be obtained by specifying how initial conditions are set before evolving towards a glass. Assuming the common choice of zero initial velocities, we have and then

 D=1a1/4 [cos(√234lna)+1√23sin(√234lna)], (39)

where we have switched to the scale factor as a time variable (recall ). Therefore, both eigenmodes contribute an overall decaying factor of with oscillations of opposite phase. The top panel in Fig. 14 shows Eq. (39) as a function of (solid), showing that the decay of perturbations due to repulsive gravity is very slow. These damped oscillations explain the behavior seen in Fig. A2 in Baugh et al. (1995).

Using Eq. (15) with the propagator in Eq. (38) and Gaussian Poisson initial conditions we can derive the behavior of the third moment as evolution proceeds towards a glass,

 ⟨δ3⟩glass = −3D2299√a[52+(23a1/4−75)cos(√234lna) (40) +√23(17a1/4−9)sin(√234lna)] σ4P,

where denotes the initial Poisson variance. This result is shown as a function of in the bottom panel in Fig. 14 (solid). Note that (mostly negative) skewness is developed as evolution to the glass takes place, therefore although the amplitude of perturbations mostly decreases as things evolve to a glass, skewness is generated. Unless one evolves by many () scale factors this can potentially enhance the problem of transients. The reason is that if the glass itself has negative skewness, this adds to the incorrect skewness that initial conditions may have, as we now discuss. This issue does not arise for grid initial conditions.

Particles are displaced from their glass positions by the displacement vector so that their positions are given by Eq. (29), and the induced density perturbations then obey

 (1+δ(q))d3q=(1+δ(x))d3x, (41)

where are the density perturbations due to the distribution of starting positions. These vanish for a grid below the Nyquist frequency, but are non-zero for a glass. Then, linearizing, the total perturbation imposed on a set of particles is

 δ(x)=(1J−1)+δ(q)≡δic(x)+δglass(x), (42)

where the first term is the standard initial condition perturbation generated by the displacement field, with given by Eq. (32), and the second term is the glass perturbation, the outcome of the process described by Eqs. (38-40). Since the two perturbations are independent, requiring that the glass skewness is negligible (and thus has no impact on transients) means to require it to be small compared to the Poisson noise of the particles (which typically is smaller than the physical perturbations imposed by the initial condtions). Then one must simply require

 ⟨δ3⟩glass≪σ4P, (43)

which according to Fig. 14 requires in the evolution towards a glass.

In the desire for faster convergence towards a glass, the Hubble drag term in the equations of motion is sometimes increased to avoid oscillations (Couchman et al., 1995). Making such a modification in Eq. (37) to allow for a generic drag of magnitude gives,

 Ωab≡[0−13/2ν], (44)

which leads to propagator eigenmodes time dependence as with , and thus the requirement of no oscillations in approaching a glass gives . Taking as an example, we obtain the growth factor and third moment as before, shown in Fig. 14 as dashed lines. We see from this that although the convergence to a glass is improved in the growth factor sense, the generation of skewness is more significant, and therefore one must wait roughly the same as with to have negligible impact.

Summarizing, generating a glass with an overall expansion factor by should be a safe configuration of points for starting cosmological simulations.
{"extraction_info": {"found_math": false, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.9570750594139099, "perplexity": 1379.6989434651132}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2021-31/segments/1627046153392.43/warc/CC-MAIN-20210727135323-20210727165323-00339.warc.gz"}