Introduction

Cosmological N-body simulations are a well-established technique for understanding nonlinear structure formation. A common technique for N-body simulations is model a single fluid, corresponding to a combination of cold dark matter (CDM) and baryons, under gravity, under the approximation that these two components trace each other. Hydrodynamic simulations follow baryons using a separate particle species, but frequently use the same initial transfer function for each species [e.g. 1].

However, at early times the baryons couple to radiation, while the CDM does not. This both induces the baryon acoustic oscillation (BAO) peak and reduces the clustering of baryonic matter on scales less than the comoving horizon scale at recombination. Structure formation reduces this relative difference between baryons and CDM, as each species couples equally to the gravitational potential. By z = 0 the power spectra of baryons and CDM differ by less than 1%. Nevertheless, at higher redshift they differ substantially, by up to 8% at z = 10 and 3% at z = 2 for k > 0.02 k/Mpc.

As observational probes of the Universe reach ever higher redshift, these differences may become more relevant. However, naively modelling separate transfer functions for each species in simulations does not reproduce a key prediction of linear theory, the evolution of the offset in power between the CDM and baryon species [2,3]. We show that using two offset grids for the baryons and CDM causes the relative differences in power to be over-estimated due to a spurious growing mode. Even with separate species, simulations correctly evolve the total matter spectrum and thus may be used to produce, for example, accurate galaxy halo catalogues. However, some cosmological probes, such as the Lyman-α forest [4], 21-cm power spectrum [5], the earliest structures [6][7][8] or even globular clusters [9] are sensitive specifically to gas at high redshift. An accurate cosmological analysis of these probes should correctly reproduce the linear theory evolution of each species.

In this paper we discuss several methods of achieving this. Our preferred resolution is to initialize the (pre-displacement) baryons using a Lagrangian glass rather than a grid offset from the CDM. A Lagrangian glass initializes particles using a homogeneous quasi-random tessellation of the mass density derived by minimizing the gravitational potential [10]. Other solutions are over-sampling the CDM relative to baryons by a factor of Ω CDM /Ω b (so that particle of each type has the same mass), or using an adaptive gravitational softening for gas particles. We evaluate and compare each method. We also discuss how linear perturbation theory can partially clarify the discrepancy that occurs with two offset grids. Finally, we show that using species-dependent transfer functions affects the Lyman-α forest at the 5 -10% level on scales measured by the BOSS survey [4].

Several earlier works have treated aspects of the same question. Refs [2, 3] showed that linear theory can be reproduced using adaptive softening lengths. With adaptive softening lengths, the gravitational softening of the baryon particles is set proportional to their SPH smoothing length. This effectively disables the short range gravitational force at high redshift, when the Universe is homogeneous and SPH smoothing lengths are of order the mean interparticle spacing.

This procedure is sometimes justified on the basis of the Jeans length of a pressuresupported cloud [11] We emphasize that on cosmological scales this justification does not apply as pressure forces are negligible on these scales. The specific adaptive softening length chosen, the SPH smoothing length, thus cannot be justified on the basis of hydrodynamics. The coupling effect is purely gravitational, and persists even if (as in our simulations and those of Ref. [3]) hydrodynamical pressure forces are disabled. In practice therefore, adaptive softening has an effect similar to uniformly increasing the softening only at high redshifts where the problem is most acute. Since the small-scale gravitational force is disabled at high redshift, adaptive softening strongly suppresses the small-scale clustering of the baryons in low density regions outside galaxies. This may make it unsuitable for some applications, in particular simulations of the Lyman-α forest. Nevertheless, adaptive softening is widely used in recent simulations [e. g 12]. 1An alternative resolution to this problem was proposed by Ref. [14]. They showed that the desired relative power between baryons and cold dark matter could be achieved by initializing the particles using two independent Lagrangian glasses, rather than the traditional offset grids. Several of our results are implicit in Ref. [14], and in related work discussing the effect of the small-scale supersonic baryon streaming velocity [15]. However, they did not consider large scales and thus did not include the effect of species and scale dependent velocity transfer functions as we do here. In particular, on the scales we consider the relative velocity baryon-CDM streaming velocity [16] is naturally included as part of the initial velocity transfer function.

Other possibilities for the initial particle distribution have also been proposed. Ref. [17] used a quaquaversal tiling of three dimensional space, which we do not consider as it restricts the total number of particles to a power of two. Ref. [18] used a capacity constrained Voronoi tessellation, which preserves most of the good properties of a Lagrangian glass. We do not use it here only because it is not obvious how to generalize it to two fluids. Both of the above are examples of other quasi-random, homogeneous tessellations of space, like a glass file. Finally, Ref. [19] avoid the problem (but only at z = 0) by "backscaling", where each species is scaled separately by the requisite factor so that the final simulation output exactly matches the linear theory prediction.

The rest of this paper proceeds as follows: in Section 2 we discuss in detail our techniques for initialising cosmological simulations. In Section 2.1 we review the formation of a glass file. Section 2.2 explains how we set the initial displacements and velocities of the particles. Section 2.3 describes our example simulation suite. Our results for the matter power spectrum are described in Section 3. Section 3.1 reproduces the existing discrepancy with linear theory. Section 3.2 shows that there is good agreement with linear theory for our preferred strategy of a Lagrangian glass for baryonic particles. Section 3.3 examines other successful strategies, comparing their strengths and weaknesses. Section 4 discusses and explains our results, noting the mechanism by which each successful simulation strategy works. Section 5 shows the effect on the Lyman-α forest of a species-dependent transfer function. Finally, Section 6 concludes.

Methods

In this Section we describe our methods for initializing cosmological simulations, as well as the simulations we run. The basic algorithm of our preferred setup is similar to that presented in Ref [14] and is as follows:

1. Produce a (homogeneous) distribution of CDM particles arranged in a regular grid.

2. Produce a (homogeneous) distribution of gas particles from a Lagrangian glass.

3. Evolve the mixed distribution of both types of particles using a reversed gravitational force, as when generating a Lagrangian glass. This avoids random close juxtapositions of CDM and baryon particles.

4. Set the velocity of the baryon and CDM particles using the scale-dependent velocity transfer function specific to each species.

5. Displace the baryon and CDM particles using the displacement transfer function specific to each species.

Section 2.1 describes Steps 1-3, which generate a homogeneous distribution of particles with minimal residual power proportional to k 4 . These are used as our pre-displacement initial sampling of matter. In Steps 4-5, this homogeneous particle distribution is displaced so as to be a Gaussian random field with the velocity and position power spectrum given by cosmological perturbation theory. These steps are described in Section 2.2.

Grid and Glass Particle Distribution Generation

We generate glass particle distributions following Ref. [10]. A gravitationally interacting N-particle system is evolved with the sign of gravity reversed, so that each particle feels a repulsive gravitational force from every other particle. We have implemented this gravitational force evolution directly into our initial conditions code and it thus uses the particle-mesh algorithm from our N-body code, MP-Gadget. The initial particle distribution before glass generation is a regular grid with the positions of each particle randomized over three mean inter-particle spacings. A symplectic leap-frog is run with a kick-drift-force-kick arrangement for 14 timestepsfoot_1 , the smallest number of steps that produced a small-scale noise power spectrum in good agreement with the expected k 4 shape. The velocity kick includes a damping term proportional to the velocity to avoid particle positions oscillating around the minimum. The velocity produced by the glass evolution is zeroed after the glass generation is complete.

This procedure can be applied separately to both CDM and baryons, generating two uniform homogeneous gas distributions. Although the distribution of each individual particle species is homogeneous, there is no guarantee that this property also applies to the total particle distribution. The original benefit of a glass, low residual force between all particles, does not hold for two combined glasses [14] as particles of one species may occasionally, by chance, be initialized extremely close to a particle of a different species. This causes two types of problems: the first is that close pairs of particles causes the computation of the gravitational force to be computationally expensive. The second is that there will be an initial over-density around the particle pair, which may source numerical structure growth.

We investigated generating a single coherent glass file containing both species of particle. In this glass file both particle species are placed at random and then evolved with reversed gravity to a single, coherent, homogeneous distribution. However, this contained no mechanism for keeping the distribution of each species homogeneous, and thus produced a density map which had over-densities of baryons exactly compensated by under-densities of cold dark matter. This in turn led to a very noisy power spectrum estimation as not all modes were fully sampled. Furthermore, in a realistic hydrodynamical simulation it would have caused, for example, the formation of dark matter halos containing no baryons.

Ref. [14] proposed a resolution to reduce chance alignments of two particle species. After generating two separate homogeneous glass files, they ran a single reversed-gravity timestep on the combined particle distribution. This ensures that particles which are exceptionally close by chance are disrupted, but does not disrupt the overall glass distribution of the simulation. We perform 14 reversed gravity timesteps on the combined particle distribution, so that the overall particle distribution is homogeneous. As each species is roughly homogeneous, the residual gravitational force for each species is small and so for most particles the distribution is not stongly perturbed. The number of reversed gravity steps is likely excessive, but once homogeneity is achieved more steps have no effect. We thus use 14 timesteps as a careful default for larger simulations than our tests.

Lacking the symmetry of a regular cubic grid, glass files produce spurious noise proportional to k 4 [20]. Since our use of a glass distribution is solely designed to reproduce the linear theory offset between baryons and cold dark matter, we can reduce the glass noise by using a regular grid for the cold dark matter. 3 Glass noise is still present, but at a reduced level because baryons are only ∼ 1/6 of the total matter. We also show that in practice hydrodynamical forces suppress the small scale power in baryons, reducing glass noise still further. For realistic particle loads we found that the noise from the glass quickly became subdominant to the physical power spectrum from structure formation, but we use a regular grid for the CDM nevertheless.

Initialization of Particle Velocities and Displacements

Key to reproducing linear theory is the correct initialization of particle velocities and displacements. Since there are several approximations current in the literature, especially for the velocity transfer functions, we here review the linear results.

We initialize both displacements and velocities directly using the transfer functions from the linear Boltzmann code CLASS [21], corresponding to the solutions of the equations of linear cosmological perturbation theory. We apply a species-independent boost to the velocities so that the simulation takes place in a constant-time Newtonian frame. This is because the CLASS transfer functions are generated in synchronous gauge, a frame comoving with the cold dark matter, and we must convert to the Newtonian frame of the simulation code. Following the notation of Ref. [22], we set:

where ḣ and θ b are set using the relevant columns in the CLASS transfer functions. h is the synchronous gauge density perturbation, and θ b the baryon velocity dispersion. Particle displacements are set using the columns for δ CDM,b . Note that each of these quantities is fully scale-dependent, so that the scale-dependence of the velocity perturbations is automatically included without having to numerically differentiate the growth factorfoot_3 . We do not use the common Zel'dovich approximation, which is to assume a scaleindependent proportionality of the velocity transfer function to the displacement transfer functions

where x is the particle displacement and D 1 (a) is the growing mode total matter linear growth factorfoot_4 . As growth function for individual species is not the total growth function, the Zel'dovich approximation introduces an inaccuracy in the relative velocities of the baryons and CDM. We instead use the perturbation theory variables described above directly. Once the displacement field is generated it is interpolated onto the particle grid using cloud-in-cell. We do not use second-order Lagrangian perturbation theory [23] as the full species specific velocity and displacement second order transfer functions (including baryon-CDM cross-terms) do not exist in the literature. We defer this to future work.

Example Simulations

Our primary test simulations are performed with 2 × 256 3 particles and box sizes of 300 Mpc/h (comoving). Initial conditions are generated at z = 99 and we compare to linear theory power spectra and transfer functions generated using CLASS [21]. We assume a periodic box and a ΛCDM cosmology with Ω M = 0.3, h = 0.67, n s = 0.97 and σ 8 = 0.8. We have set Ω b = 0.05 so that Ω CDM /Ω b = 5. For the majority of our tests the only force included is Newtonian gravity. Even though we include baryon particles, we generally disable pressure forces. Exceptions are a single simulation to show pressure effects, and our Lymanα forest simulations. When enabled, hydrodynamics is followed using Smoothed Particle Hydrodynamics with a cubic kernel [24]. Except for our simulations of the Lyman-α forest, cooling and star formation are not included.

Our simulations are performed with the N-body code MP-Gadget and initial conditions were generated with MP-GenIC [25]. MP-Gadget is descended from Gadget-2 [24], with heavy modifications for scalability and force accuracy. In particular, as described in Ref. [26], the force tree is rebuilt completely every timestep, rather than nodes being moved according to the average particle velocity. We have also replaced the force kernel which smoothes the short-range tree force into the long-range particle-mesh force. In Gadget-2 a Gaussian error function is used. For MP-Gadget we instead estimate the smoothing kernel by comparing the particle-mesh force to a simple direct N 2 N-body force evaluation, following [27].

Our full set of simulations is shown in Table 1. To demonstrate that we reproduce the results of the literature, TWOGRID is a simulation with species-specific transfer functions and a homogeneous, pre-displacement particle distribution of two offset grids. Grids are offset equally in each Cartesian direction by factors weighted by the contribution of each particle species to the total matter density, so that the total center of mass is unchanged. We have confirmed explicitly, using a small simulation with zero initial displacement, that this is a stable configuration. The ADAPTIVE simulation used the same initial conditions but evolved them using an adaptive gravitational softening set to match the SPH smoothing length.

We performed simulations (HALFGLASS) using the setup described in Section 2.1, with pre-displacement CDM particles on a regular grid and pre-displacement baryon particles on a Lagrangian glass. To confirm that our results continue to hold on larger scales we performed a simulation (BIGGLASS) with a box of 1 Gpc/h (a scale large enough that the baryon and CDM transfer functions are approximately equal on scales of the box) and 2 × 768 3 particles. To check convergence with resolution, we performed a simulation with 300 Mpc/h and 2 × 512 3 particles (HIRESGLASS). To confirm the effect of pressure forces, HYDROGLASS uses the same parameters as HALFGLASS but enables hydrodynamics for the baryonic particles.

The UNDERSAMP simulation undersamples the baryon particles so that each CDM and baryon N-body particle has the same mass. For this simulation the pre-displacement CDM and baryon particles were initialized on offset grids, as for the TWOGRID simulation, but the mean inter-particle spacing of the baryon particles was 5 times wider than the CDM.

Finally, we performed two simulations adapted for the Lyman-α forest, in a 120 Mpc/h box with 2 × 512 3 particles. These parameters were chosen so that the box was large enough to still include some portion of the BAO, while having sufficient resolution to roughly model the forest. One simulation (FOREST) used species specific transfer functions. A second (TOTFOREST) contains CDM and baryons which are both initialised using the total matter transfer function. Both simulations include radiative cooling following [28], a self-shielding correction from [29] and a simplified star formation criterion which immediately turns all gas with an over-density ∆ > 1000 and a temperature T < 10 5 K into stars [30]. A uniform meta-galactic ultra-violet background was assumed following Ref. [31].

Results

In this Section we describe the results of our simulations. In Section 3.1, we discuss the discrepancy with linear theory. In Section 3.2 we show that our preferred setup of a Lagrangian glass for baryonic particles produces good agreement with linear theory. Section 3.3 compares the various simulation strategies which successfully reproduce linear theory.

Offset Grid Simulations

Figure 1 shows baryon and cold dark matter power spectra for the TWOGRID simulation, compared to the output of CLASS at the same redshift. This simulation was initialized using the traditional method of two offset grids of particles. Figure 1 illustrates the problem we aim to solve. Although the total matter power spectrum obeys the linear predictions (right panel), the DM power spectrum is over-predicted and the baryon power spectrum under-predicted, leading to an error in the relative power spectrum. 6 Disabling the short-range force entirely does allow the simulation to reproduce linear theory [3], so it is tempting to attribute the discrepancy to a simple force inaccuracy in the Gadget short-range force. We have confirmed that our results are independent of the force accuracy of the simulation. We performed a simulation where the short range force was evaluated using an N 2 pairwise gravity solver a simulation where the tree opening angle was changed and a simulation where the split scale between the short and long-range forces was increased. Although the total power in the box in some cases changed on small scales (as expected given the force accuracy), the power ratio between the baryons and cold dark matter was changed by < 0.5%. As a further check, we have compared the gravity force from MP-Gadget to an independently written simple O(N 2 ) particle solver and found good agreement.

We checked that our results were unchanged when all particles had the same timestep and when the long-range PM timestep was 4 times shorter, ruling out an effect due to adaptive timestepping or the length of the timesteps. We also checked that the results of Figure 1 persist when the box size and particle number were increased. Finally, we have performed the same simulation using the Gadget-2 code, with consistent results, see Appendix A. Some of the scales simulated are already non-linear by z = 2; Figure 1 shows a non-linear scale of k = 0.25 h/Mpc. There therefore exists the (unlikely) possibility that the discrepancy is a physical effect resulting from non-linear growth. To rule this out, we performed a simulation with suppressed initial power, so that σ 8 (z = 0) = 0.008. This simulation gave a similar result to the TWOGRID simulation, but with increased noise on small scales.

Our broad results, that two simulation grids lead to inaccurate power on large scales and the inaccuracy depends on the softening length, match those of Ref.

[3] (AHA13). However, some aspects of our simulations differ. In particular, the simulations of AHA13 produce more power in the baryon component than linear theory predicts. For a similar simulation, we find less power than predicted by linear theory. 7 However, the size of the discrepancy found by AHA13 depends strongly on redshift (Angulo, private communication). We therefore suspect some other issue was at play in the earlier work, possibly an implicit choice of a different gauge. Although the physical evolution must be the same on sub-horizon scales for all well-specified gauges, numerical artifacts need not be. Note that the disagreement is 6 The relative power drops at k > 1 h/Mpc as k approaches the scale of the particle grid at 5 h/Mpc. Harmonics of the grid scale begin to be important at the percent level by k ∼ 1 h/Mpc. 7 Figure 1 of AHA13 also shows that the magnitude of the discrepancy changes non-monotonically with the smoothing length. For example, the largest P b /P cdm occurs with a PM-only code and a smoothing length of 300 kpc/h, similar to the mean inter-particle spacing of 500 kpc/h, whereas AHA13 show that a simulation with a softening length of 100 kpc/h is much closer to linear theory. We performed a similar PM-only simulation using a PM grid cell size of half the mean interparticle spacing and successfully matched linear theory. We have never observed P b /P cdm > 1 as shown in that figure. We note that if their PM2 and Tree2 simulations were accidentally showing P cdm /P b then the discrepancy would be monotonic with smoothing length and the PM2 simulation would nearly agree with linear theory and our results. The simulations shown are no longer available, so this is impossible to confirm (Angulo, private communication).

Section 2.3. We generate 32, 000 Lyman-α absorption spectra using Ref. [33]. Sightlines are placed at random positions in the simulation box. The flux power spectrum shown is the 1D power spectrum of the flux, computed along the line of sight to the quasar.

Figure 6 shows our results. The left panel shows the ratio in the flux power spectrum between the simulation using species specific initial transfer functions and the simulation using the total matter transfer function. As might be expected, the reduction in the baryonic power spectrum at high redshift leads to a reduction in the flux power spectrum. The effect is larger than that on the matter power spectrum. This may be because the Lyman-α forest at high redshift is sensitive to gas at over-density 10 -100, at scales where the dark matter has only recently started to grow non-linearly. Since the power in gas is suppressed, it will start non-linear growth later and the difference between the two species will be increased (as non-linear growth is faster) 10 . Ultimately virialization erases differences between gas and dark matter completely.

At z = 2.2 the power in the Lyman-α forest is actually enhanced. This may be because the Lyman-α forest is sensitive to redshift space effects. In order for the gas potential to catch up to the dark matter potential, there should be a higher gas velocity power spectrum. Redshift space effects are more important at lower redshifts because the forest probes slightly larger overdensities.

Realistic forest analyses marginalize over the uncertain thermal history of the intergalactic gas as well as the mean flux of the Lyman-α forest power spectrum [e.g. 4]. The mean flux controls the overall amplitude of the flux power spectrum. The different density distributions in the TOTFOREST and FOREST simulations produce subtly different recombination rates and neutral fractions, and thus have moderately different thermal histories and mean fluxes. In order to establish how much this effect is driving the differences in the flux power spectrum observed in Figure 6, we rescaled each spectrum by a constant factor so that they have the same mean transmitted optical depth, τ = 0.0023(1 + z) 3.65 [35]. We also fit a power law temperature-density relation and rescaled the temperature of every particle in the box by a density-dependent factor (∼ 2%) so that the slope and intercept of the temperature-density relation were the same in both simulations. The ionization equilibrium neutral fraction was recomputed using the new particle temperatures and new spectra generated.

The right panel of Figure 6 shows the flux power spectrum ratios after this rescaling. The effects are now much smaller, on the order of 5%. While comparable to the 1σ statistical errors on the BOSS DR9 1D flux power spectrum, this is larger than the 2 -4% statistical error achieved in the newest DR14 data release of Ref. [36].

We checked the effect of mass resolution on Figure 6 using two simulations with 60 h/Mpc boxes and 2 × 512 3 particles. The raw flux power spectra ratios shown in the left panel increased by a uniform 5%, reflecting a change in the relative mean flux. Once the spectra were rescaled to achieve the same mean flux, the difference was much smaller and consistent with noise for z >= 3. At z = 2.2 the higher resolution simulation exhibited a reduction in power by ∼ 1%. For the forest, the effect of finite resolution should be strongest at high redshift [37] and the effect of finite box size strongest at low redshift. We therefore believe that this small change is largely driven by the smaller box and thus Figure 6 shows only the lower-resolution, larger box size, simulation.

Conclusions

We have shown that the inability of most cosmological N-body simulations to reproduce the linear theory prediction for the offset between baryons and CDM arises from their use of two offset grids for baryon and CDM particles. We identify a spurious growing mode in this setup proportional to a 1/2 which dominates for z ≤ 9.

We demonstrate several possible mechanisms for resolving this issue. Our preferred solution is to use a Lagrangian glass for the baryon particles. Chance close pairs of CDM and baryon particles are mitigated by reversed gravity timesteps run on the total particle distribution. We prefer this solution as it fully resolves the discrepancy with linear theory without compromising the simulation resolution. Other considered mechanisms are adaptive gravitational softenings and under-sampling baryon particles. We show that both reduce the effective resolution of the simulation with fixed particle number. Adaptive gravitational softening in particular means the effective force resolution is reduced to the mean interparticle spacing.

Finally, we examine whether this linear theory prediction affects a specific observable, the Lyman-α forest 1D flux power spectrum. The Lyman-α forest is a natural place to look for observational consequences because it is sensitive to neutral hydrogen gas and measures high redshifts. We find that it does at the 5 -10% level. This is larger than the 2 -4% statistical error achieved in the current DR14 BOSS data by Ref. [36], and thus has the potential to affect cosmological parameter constraints from that dataset [e. g 38]. Given the complexity of the Lyman-α likelihood function, the exact parameters affected and the magnitude of the shift is not clear, and we defer this question to future work using a larger simulation suite.

We have not investigated the impact on small-scale constraints on the warm dark matter mass [39] as our current simulations do not resolve the flux power spectrum on scales of k = 0.01-0.1 s/km (k = 0.6-6 h/Mpc at z = 4). The statistical errors on these scales are still of order 10%, so it is likely that the systematic we discuss here has minimal effect. However, our work demonstrates that a species dependent transfer function could be important for future constraints. There may also be observational implications for future high redshift 21cm or intensity mapping constraints.