FAST GENERATION OF ENSEMBLES OF COSMOLOGICAL N-BODY SIMULATIONS VIA MODE RESAMPLING

Altering the amplitudes and phases of the large-scale modes of the mass density in a cosmological volume has two dominant effects:

• 1.  
  sub-volumes are expanded or compressed, and
• 2.  
  sub-volumes gravitationally evolve with an effective matter background density that is enhanced or diminished,

depending on where the large-scale modes become more or less overdense, respectively (Frenk et al. 1988; Little et al. 1991; Tormen & Bertschinger 1996; Cole 1997; Angulo & White 2010). We follow Cole (1997) to model both of these effects. To model the expansion and compression of sub-volumes, we perturb the particle positions using the Zel'dovich displacements (Zel'dovich 1970) calculated from the large-scale modes, δ_L,

$$\begin{equation} \mathbf {d}_{L}(\mathbf {x},a) \equiv \int \frac{d^3k}{(2\pi)^3} \frac{\delta _L(\vec{k},a)\vec{k}}{k^2}\, e^{i\vec{k}\cdot \mathbf {x}}. \end{equation} \tag{ 1 }$$

The expansion and compression of sub-volumes is analogous to the frequency modulation of a radio wave. The large-scale modes to be added are a new signal that modulates the sizes of sub-volumes in the matter density field; expressed to first order by the Zel'dovich displacements in Equation (1).

Compressed (expanded) regions also evolve with an increased (decreased) effective matter density, altering the growth rate of structure. Because the growth rate increases monotonically with time, we model the changing effective local matter density by finding new times $a^{\prime }(\mathbf {x})$ at which the linear growth rate matches the growth rate in a universe with matter density $\Omega _{m}\left(1+\delta _{L}(\mathbf {x}, a)\right)$ (with Ω_Λ kept fixed). The perturbed scale factors $a^{\prime }(\mathbf {x})$ are then defined by the relation

$$\begin{equation} D(a^{\prime },\Omega _m) \approx D(a, \Omega _m \left(1+\delta _L(\mathbf {x},a)\right)), \end{equation} \tag{ 2 }$$

where D(a, Ω_m) is the linear growth function. In practice, we solve Equation (2) numerically to get a' as demonstrated in Figure 1.

Zoom In Zoom Out Reset image size

Combining the two effects of expansion/compression and time perturbations, the model for adding large-scale fluctuations to the mass density field is

$$\begin{equation} \delta (\mathbf {x}, a) = \delta _{L}(\mathbf {x},a) + \delta _{S}(\mathbf {x}^{\prime }, a^{\prime }), \end{equation} \tag{ 3 }$$

where $\delta (\mathbf {x},a)\equiv \rho (\mathbf {x},a)/\bar{\rho }- 1$ is the mass overdensity that would be measured in the all-modes simulation, $\delta _{S}(\mathbf {x},a)$ is the mass overdensity measured in the small-modes simulation, δ_L has zero power for $|\vec{k}| > k_{\rm thresh}$, a' is defined by Equation (2), and $\mathbf {x}^{\prime } \equiv \mathbf {x}+ \mathbf {d}_{L}(\mathbf {x},a)$. Note that for this definition of $\mathbf {x}^{\prime }$ the Zel'dovich displacements are added to the particle positions prior to application of the time perturbation (i.e., at the positions at time a, not at a'). We will discuss this choice below.

When |a' − a| ≪ 1 (which requires $\langle |\delta _L(\mathbf {x})|^{2}\rangle ^{1/2} \ll 1$) the first-order correction to δ_S in a' should be sufficient for adding δ_L to the density field. Expanding both sides of Equation (2) to first order in (a' − a) and δ_L,

$$\begin{eqnarray} D(a^{\prime }) &\approx D(a) + (a^{\prime }-a)\frac{{\it dD}}{da}\nonumber \\ D(a,\Omega _m(1+\delta _L(\mathbf {x},a))) &\approx D(a,\Omega _m) + \Omega _m\delta _L(\mathbf {x},a)\frac{{\it dD}}{d\Omega _m}.\quad \end{eqnarray} \tag{ 4 }$$

Solving for a' gives

$$\begin{equation} \frac{a^{\prime }}{a} -1\approx \delta _L(\mathbf {x},a) \left(\frac{d\ln D}{d\ln \Omega _{m}} \left/ \frac{d\ln D}{d\ln a}\right.\right), \end{equation} \tag{ 5 }$$

which shows that the perturbed scale factor is roughly proportional to the large-scale density perturbations being added. Using Equation (5) in the Taylor expansion of Equation (3) gives

$$\begin{eqnarray} &&\delta (\mathbf {x},a) \approx \delta _{L}(\mathbf {x},a) + \delta _{S}(\mathbf {x}^{\prime },a)\nonumber\\ &&\quad +\,\delta _{L}(\mathbf {x},a)\, \frac{\partial \delta _{S}(\mathbf {x}^{\prime },a)}{\partial \ln a}\, \left(\frac{d\ln D}{d\ln \Omega } \left/ \frac{d\ln D}{d\ln a}\right.\right). \end{eqnarray} \tag{ 6 }$$

This suggests that the addition of large-scale modes can be quickly computed given an estimate of the nonlinear growth rate of the small-modes simulation. With a single snapshot, the growth rate can be estimated from the continuity equation,⁴

$$\begin{equation} H(a)\frac{\partial \delta (\mathbf {x},a)}{\partial \ln a} = -\frac{1}{a}\nabla \times \left((1+\delta)\mathbf {v}\right), \end{equation} \tag{ 7 }$$

and, using Equation (3),

$$\begin{equation} \frac{\partial \delta _S(\mathbf {x},a)}{\partial \ln a} = \frac{\partial \delta (\mathbf {x},a)}{\partial \ln a} - \frac{d\ln D}{d\ln a}\delta _L(\mathbf {x},a). \end{equation} \tag{ 8 }$$

So, given particle positions and velocities at scale factor a for a simulation that has no large-scale modes, Equation (6) provides a method to add new large-scale mode realizations when the large-scale density perturbations are small.

According to Equation (5) the range of a' − a is proportional to the range of $\delta _{L}(\mathbf {x}, a)$. Because $\delta _{L}(\mathbf {x},a)$ is Gaussian distributed by construction, we can determine the range of δ_L from its variance,

$$\begin{equation} \sigma ^{2}\left(\delta _{L}(\mathbf {x},a)\right) \equiv \int _{k_{F}}^{k_{\rm thresh}} d\ln k\, \Delta ^{2}(k)\, W_{{\rm top\hbox{-}hat}}^{2}(k R), \end{equation} \tag{ 9 }$$

where k_F ≡ 2π/L_box is the fundamental frequency of the simulation box, Δ²(k) ≡ k³P(k)/(2π²) is the theoretical dimensionless matter power spectrum, and W_top-hat(kR) is the Fourier transform of a spherical top-hat window of radius R. The scale R is roughly given by either the pixel scale of the gridded density field or the mean interparticle spacing, whichever is larger. For the simulations we use L_box = 1000 h⁻¹ Mpc, k_thresh = 8k_F ≈ 0.05 h⁻¹ Mpc, and, with 256³ particles, R ∼ 2π/(256k_F). Evaluating Equation (9) gives σ ≈ 0.17, which is still less than 0.2 as we required when choosing k_thresh. But, as shown in Figure 1, the variations in $\delta _{L}(\mathbf {x},a)$ greater than about 2σ lead to |a' − a| > 1 making Equation (6) insufficiently accurate. In our test simulations δ_L actually varies over ±4.3σ ≈ ±0.87, giving a range of perturbed scale factor 0.26 ≲ a' ≲ 4.4 (shown by the dotted lines in Figure 1).

In Figure 2 we show the maximum value of |a' − a| as a function of the simulation box size while keeping k_thresh/k_F fixed at 8 and, as in the preceding paragraph, assuming maximal variations in $\delta _{L}(\mathbf {x},a)$ of 4.3σ. For L_box ∼ 3000 h⁻¹ Mpc, as in the recently completed MXXL simulation,⁵ the variation in a' − a is less than 0.2 and the approximate formula in Equation (6) could be sufficient for adding the large-scale modes to the small-modes simulation. To further assess the validity of Equation (6), we show a histogram of the values of the third term in Equation (6), normalized by $\delta _S(\mathbf {x}^{\prime },a)$, on a 512³ grid in our 1 h⁻¹ Gpc small-modes simulation in Figure 3. The distribution of grid cell values is strongly peaked around zero, but has wide tails extending to values of several hundred. This again illustrates that Equation (6) is a bad approximation for a (1 h⁻¹ Gpc)³ volume. However, if we scale the bin values on the lower horizontal axis by the ratio of the maximum perturbed scale factor, max(a' − a), in a (3 h⁻¹ Gpc)³ simulation to the maximum (a' − a) in our simulation, we get the axis on the top of Figure 3. For a (3 h⁻¹ Gpc)³ simulation, 92% of the grid cells of the third term in Equation (6) have values less than one. So, Equation (6) may be a good enough approximation for adding large-scale modes to the MXXL or other simulations with similarly large volumes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For box sizes smaller than several Gpc on a side, we need an algorithm that accurately implements the growth matching condition of Equation (2) to compute the density $\delta _S(\mathbf {x}^{\prime },a^{\prime })$ in Equation (3). For computational expediency and for later comparison with Takahashi et al. (2009), we have limited our test simulations to L_box = 1000 h⁻¹ Mpc. Therefore, we have not tested the accuracy of Equation (6). Instead, because our test simulation parameters yield |a' − a| > 1, we apply the model of Equation (3) by adjusting the positions of each particle in the simulation. If $\mathbf {x}_{i}(a)$ is the position of the ith particle at scale-factor a in the small-modes simulation, the model from Equation (3) with modes δ_L added to the simulation can be rewritten for each particle position as

$$\begin{equation} \mathbf {x}_{i}^{\prime }(a^{\prime }; \delta _{L}) \equiv [ \mathbf {x}_{i}(a) + \mathbf {d}_{L}(\mathbf {x}_{i}(a), a^{*})]_{a = a^{\prime }}, \end{equation} \tag{ 10 }$$

where a* is the scale factor of the target time where the full density field (with new large-scale modes) is desired. That is, we first apply the same Zel'dovich displacement to the particle positions in all the snapshots (evaluated at the unperturbed positions in each snapshot). Then, we interpolate between these modified snapshots to find the particle position at perturbed time $a_{i}^{\prime } \equiv a^{\prime }(\mathbf {x}_{i} + \mathbf {d}(\delta _{L}(\mathbf {x}_{i}(a)), a^{*}))$. To find the particle positions at all values of a' we must save simulation snapshots at times covering the range of a' values and spaced sufficiently closely in a so that the interpolation of positions is accurate. This requires running the simulation far into the future (a > 4 for our test simulations) and saving of order tens of snapshots, which creates both an extra computational and also a disk storage burden over the approximation in Equation (6).

The steps of our full algorithm for adding large-scale modes to the small-scale density are then

• 1.  
  choose a threshold scale k_thresh defining the boundary between large- and small-scale modes and a target time a* where the resampled density is desired;
• 2.  
  calculate the range of perturbed scale factor a' defined by Equation (2) about a* where the simulation snapshots will be required;
• 3.  
  run an N-body simulation in a periodic volume with nonzero Fourier modes in the initial conditions only for $|\vec{k}| > k_{\rm thresh}$ and save snapshots over the required range of a' that are spaced sufficiently close together to allow for accurate interpolation of particle positions;
• 4.  
  calculate a new realization of large-scale modes δ_L(a*) on a grid;
• 5.  
  calculate the Zel'dovich displacements defined in Equation (1) from δ_L;
• 6.  
  apply the Zel'dovich displacements to each snapshot of the small-modes N-body simulation; and
• 7.  
  move each particle in the small-modes simulation to its location at time $a^{\prime }(\mathbf {x}(a)+\mathbf {d}_L)$ as defined by Equation (2), interpolating between snapshots when necessary.

The resulting particle positions give a close approximation to the particle positions that would be obtained by running an N-body simulation with the Fourier modes in δ_L present in the initial conditions (with the amplitudes suitably scaled according to the linear growth function).

The two operations of adding Zel'dovich displacements and perturbing the time where the particle positions are evaluated do not commute when the density field has evolved nonlinearly from the initial conditions. To decide on how to order these operations we imagine an approximation to an N-body simulation where the initial conditions are first evolved by adding Zel'dovich displacements to the particles and then the density is scaled by the linear growth function. Reversing these operations, with the Zel'dovich move applied last, would be less analogous to the way the particle positions actually evolve in a simulation. In this latter scenario, the particles would first gravitationally collapse at rates determined by the local mean matter density and then at late times would undergo displacements to expand or compress regions in under- or overdense volumes. There is nothing physical about this ordering of operations. We have confirmed with our test simulations that the mode resampling does not work well unless the Zel'dovich move is the first step of the mode-resampling algorithm.

First we define the estimator for the covariance matrix that we use for our comparisons.

Let $\delta _{\vec{k}_i}$ denote the amplitude at grid point $\vec{k}_i$ of the discrete Fourier transform of the density field from an N-body simulation measured on a three-dimensional uniform periodic grid. Then, the standard estimator for the power spectrum of the mode amplitudes is obtained by averaging $|\delta _{\vec{k}_i}|^2$ over a spherical shell of constant radius $|\vec{k}_i|$,

$$\begin{equation} \hat{P}_q = \frac{V}{N_{k_q}} \sum _{\vec{k}_i, |\vec{k}_i|=k_q} \big|\delta _{\vec{k}_i}\big|^2, \end{equation} \tag{ 12 }$$

where $N_{k_i}$ is the number of modes that fit in the shell. We follow Takahashi et al. (2009) and set this shell thickness to 0.01 h Mpc⁻¹. If we generate N_r realizations of the density field (e.g., by running N_r identical N-body simulations except with different random number seeds to generate the initial conditions), the mean power spectrum estimator for all the realizations is

$$\begin{equation} \bar{P}_q = \frac{1}{N_r}\sum _{a=1}^{N_r}\,\hat{P}_q^{a}. \end{equation} \tag{ 13 }$$

The estimator for the covariance of the power spectrum from all the realizations is then

$$\begin{equation} \hat{C}_{ij} = \frac{1}{N_r}\sum _{a=1}^{N_r} (\hat{P}_i - \bar{P}_i)\, (\hat{P}_j - \bar{P}_j). \end{equation} \tag{ 14 }$$

Meiksin & White (1999), Norberg et al. (2009), and Takahashi et al. (2009) found that at least several hundred realizations of the matter density field are required to accurately estimate the sample covariance matrix. This means that hundreds (or thousands in some cases) of N-body simulations must be run with different random number seeds in order to obtain a power spectrum covariance matrix estimate that is sufficiently accurate for cosmological parameter estimation.

Because running hundreds or thousands of N-body simulations is often prohibitively expensive, previous studies have attempted to estimate the power spectrum covariance from a single N-body simulation by applying window functions to the gridded density field to sub-sample the simulation volume. Hamilton et al. (2006) found that because the windows are convolved with the Fourier modes of the periodic simulation box, the resulting power spectrum covariance matrix estimates are significantly biased with respect to the result obtained from an ensemble of periodic simulations. Furthermore, this bias cannot simply be accounted for by deconvolving the window functions from the power spectrum estimates, but is a result of higher order connected correlations entering the covariance matrix estimator, which Hamilton et al. (2006) call "beat-coupling."

Although direct sub-sampling of the simulation volume has been shown not to reproduce the covariance from an ensemble of simulations, the idea that different regions of a large simulation volume provide statistically independent information about the small-scale covariance is sound if we accept the ergodic hypothesis and assume our simulation is big enough. To avoid the convolution in Fourier space imposed by the configuration-space windows of Hamilton et al. (2006), we attempted to sub-sample the Fourier modes that enter the shell-averaged power spectrum estimator in Equation (13). For shells with large wavenumbers there are thousands of modes used to estimate the power spectrum. The nonlinear power spectrum covariance depends on the connected four-point correlations of the Fourier modes, but with thousands of modes on a uniform grid there must be many redundant four-point configurations that could yield independent samples of the (shell-averaged) power spectrum covariance (after appropriate scaling by the number of modes used to estimate the mean power spectrum).

The enumeration and counting of the four-point configurations in each shell is complicated, so we instead divided the surface of each wavenumber shell into uniform area HEALPix (Górski et al. 2005) pixels and then selected only one wavevector within each pixel to obtain a power spectrum estimate. For a given number of pixels covering the spherical shell, the minimum multiplicity of wavevectors in a pixel determines the number of pseudo-independent power spectrum estimates we can derive with this method. That is, if all pixels have at least n wavevectors within their boundaries, then we can generate at most n pseudo-independent realizations of the power spectrum estimates. We defined smaller (and therefore more) pixels for shells of increasing wavevector magnitude. This ensured that the multiplicity of each pixel was greater than one for a wide range in wavevector magnitudes.

We applied this sub-sampling to a small-modes simulation with 500 large-scale mode resamplings. The resulting power spectrum estimates and correlation coefficients of the estimated power spectrum covariance are shown in Figure 9. While the mean power spectrum is successfully recovered with this algorithm, the power spectrum covariance estimator has nearly zero off-diagonal terms, in stark contrast to the result from Takahashi et al. (2009; shown in red, dashed lines) for wavenumbers ≳ 0.1 h Mpc⁻¹. One explanation for this could be that our method of sub-sampling modes in the wavevector shells is not actually selecting representative samples of all the four-point configurations in each shell, which would be fixed by a more precise sub-sampling procedure.

Zoom In Zoom Out Reset image size

Because pursuing this algorithm further would significantly increase the computational complexity while limiting the range of applications, we instead turned to using multiple small-modes simulations with large-scale mode resampling to further explore the covariance matrix estimates as described in the following section.

The power spectrum variance and correlation coefficients estimated from 500 resamplings of the large-scale modes in a single small-modes simulation are shown by the blue squares in Figure 10. Note that we use the same binning in wavenumber as Takahashi et al. (2009) in all the cases shown in order to make an accurate comparison. While the small-scale variance estimate is boosted beyond the Gaussian prediction, it is systematically biased low relative to the results of Takahashi et al. (2009). At intermediate wavenumbers just larger than k_thresh, the estimated variance is only ∼20% of the true value. This is not too surprising because it is at these scales where our algorithm for adding large-scale modes is expected to be most inaccurate. This is because our approximation that longer wavelength modes simply change the effective local matter density must be wrong when the longer wavelengths are nearly equal to the smaller wavelengths to be perturbed. The modes both just larger and smaller than k_thresh may contribute to coherent structures in the mass density distribution so that resampling only some of the Fourier modes in these structures would severely underestimate the true variance on these scales. Also, we should note that the modes with wavenumbers just smaller than k_thresh are likely not strictly Gaussian distributed in the all-modes simulation, so that their true distribution and coupling to small modes is not properly represented.

Zoom In Zoom Out Reset image size

The right panel in Figure 10 shows four rows of the matrix of estimated correlation coefficients of the power spectrum. It is significant that the mode resampling generates nonzero correlation coefficients between the large and small scales (blue crosses—lower panels especially), which are identically zero in the unperturbed small-modes simulation. However, the correlation coefficients are again biased low by multiplicative factors of two or more.

We have checked that the biases in the estimated power spectrum variance and correlation coefficients do not change when using more than 500 large-scale mode realizations. Instead we speculate that this bias is real and indicates that (at least for a fixed k_thresh) the mode-resampling algorithm generates realizations of the matter density from only a subset of the full nonlinear density probability distribution. This makes sense intuitively if we realize that the mode resampling as implemented here can never generate a new set of halos in the simulation box, but only moves and merges or de-merges the halos already present in the small-modes simulation. Several previous investigations have shown that the phase correlations in the Fourier modes of the density field in a dark matter simulation are dominated by the masses and positions of the most massive halos (Watts et al. 2003; Neyrinck & Szapudi 2007). It makes sense then that the power spectrum variance on small scales would also be quite sensitive to the particular realization of the most massive halos. However, with a larger simulation volume the covariance estimates might be less dependent on the particular realization of the massive halos because many different halo configurations would be present.

If our mode-resampling algorithm is simply sampling a subset of the density probability distribution, then the bias in the power spectrum variance and correlation coefficients should decrease if we use more than one simulation. To check this hypothesis we ran 19 more small-modes simulations (for a total of 20) with different random number seeds in the initial conditions and added 500 resamplings of the large-scale modes to each simulation (for a total of 10,000 realizations of the density field). The sample covariance estimates improve considerably as shown by the red diamonds in Figure 10. The biases in both the variance and the correlation coefficients are much smaller over all scales. In addition, the extremely low dip in the variance estimates at intermediate scales when using only one simulation has disappeared entirely with 20 simulations.

Finally, we show an estimate of the covariance matrix using 20 simulations without any resampling of large-scale modes with the green crosses in Figure 10. With such a small number of simulations the statistical noise dominates in the estimator. Comparing with the estimator in 500 resamplings of each of the N-body simulations shows how much our mode-addition algorithm reduces the statistical errors.

Although with 20 small-modes simulations the bias in the power spectrum covariance estimate is much smaller than with only one simulation, a small bias is still noticeable in Figure 10. In Figure 11, we attempt to quantify the improvement in the covariance matrix estimates as more small simulations are added by plotting the rms dispersion in elements of the power spectrum variance relative to the result from Takahashi et al. (2009). We show the distribution in the variance estimates for wavenumber bins with centers in 0.295 < k/(Mpc h⁻¹) < 0.395 with the boxes in the left panel of Figure 11. The horizontal lines in the centers of the boxes denote the median of the variance estimates and the boxes' tops and bottoms denote the first and third quartiles. We have selected only high-k modes to plot in Figure 11 because the variance estimate in this range has a roughly constant bias as seen in Figure 10 and we are primarily interested in the accuracy of the covariance estimates on small scales where the deviation from the Gaussian model is most severe. We show how the normalization of the fit to the covariance dispersion changes with the number of mode resamplings per N-body simulation in the right panel of Figure 11. The reduction in the number of N-body simulations needed to give a certain error in the estimate covariance flattens around 100 resamplings per simulation.

Zoom In Zoom Out Reset image size

Takahashi et al. (2009) performed a similar analysis and found that the dispersion in their sample variance estimators was well fit by σ = 2/N_r, where N_r is the number of N-body simulations (or realizations) used to compute the estimator (shown as the dashed black line in Figure 11). We find a nonlinear least-squares fit to our variance estimates of 2/(8.5N^0.8_r), which is shown by the solid red line in Figure 11. From this we conclude that the sample variance estimator using resampled large-scale modes can achieve similar accuracy as the standard sample variance estimator but with ∼8 times fewer N-body simulations. However, it is difficult to extrapolate this convergence rate far beyond the 20 simulations we have actually performed because small changes in our choice of wavenumber bins or in our least-squares fit can cause large extrapolation errors.