THE COYOTE UNIVERSE EXTENDED: PRECISION EMULATION OF THE MATTER POWER SPECTRUM

The discovery of the accelerated expansion of the universe by Riess et al. (1998) and Perlmutter et al. (1999), now confirmed by multiple observations, has had a tremendous impact on the field of observational and theoretical cosmology. Due to our near-complete ignorance regarding the origin of cosmic acceleration, the initial task is to better characterize the nature of the underlying cause, first by observational means. The primary focus in this endeavor is to decide whether the acceleration is caused by some form of dark energy, a cosmological constant being the simplest realization of this idea, or by some modification of general relativity on large cosmological scales.

Unfortunately, carrying out controlled experiments in cosmology is impossible—therefore, information relevant to the above task must be obtained by combining inferences from different kinds of observations. These observations include the measurement of the temperature anisotropy of the cosmic microwave background (CMB), the luminosity distance relation from supernova observations, and large scale structure probes, such as baryon acoustic oscillations (BAO), weak lensing, clustering of galaxies, and the abundance of galaxy clusters. Despite the impressive power of the observations, all cosmological probes are plagued with a large number of sources of systematic error. These can be broken up into two kinds: (1) those related to the observational techniques—shape measurements in weak lensing and the determination of cluster masses being two prominent examples, and (2) those related to the uncertainties in the theoretical predictions and the (related) errors in solving cosmological statistical inverse problems. An important member of this latter class is the prediction accuracy of the matter power spectrum, the central theme of this paper.

Because robust results in cosmology must rely on combining a number of observational probes, each possessing their individual strengths and weaknesses (not all of which can be fully predicted in advance), the current strategy is to proceed with several survey missions. In the spectroscopic realm, ongoing or recently completed surveys focusing on the measurement of BAO and galaxy clustering include the Sloan Digital Sky Survey (SDSS; Zehavi et al. 2011), the WiggleZ survey (Drinkwater et al. 2010) and the Baryon Oscillation Spectroscopic Survey (BOSS; Schlegel et al. 2009) component of SDSS-III; planned next-generation missions include the Mid-Scale Dark Energy Spectroscopic Instrument (MS-DESI) and SuMIRe. Imaging surveys include SDSS, the Canada–France–Hawaii Telescope Legacy Survey (CFHTLS), the Kilo Degree Survey (KiDS; de Jong et al. 2013), the recently started Dark Energy Survey (DES), and in the future, Subaru Measurement of Images and Redshifts (SuMIRe; Sugai et al. 2012), the Large Synoptic Survey Telescope (LSST; Abell et al. 2009), and Euclid (Refregier et al. 2010). The most recent constraints from weak lensing have been obtained by the CFHTLenS team (Kilbinger et al. 2013), analysis of SDSS data (Huff et al. 2011, Lin et al. 2012), and the Deep Lens Survey (DLS; Jee et al. 2013). In the X-ray the next planned mission is eROSITA (Merloni et al. 2012), which will provide a significant update to results from prior surveys, such as the 400d cluster survey (Vikhlinin et al. 2009). In the CMB, high-resolution measurements by the Atacama Cosmology telescope (ACT; Das et al. 2011) and the South Pole Telescope (SPT; van Engelen et al. 2012) have confirmed the sensitivity of the CMB to the large scale distribution of matter (CMB lensing). Significantly, first results from the Planck mission have been recently released (Ade et al. 2013). For a comprehensive review of the current state-of-the-art (pre-Planck) we refer the reader to Weinberg et al. (2013).

With the advent of ambitious large-scale observational programs, emphasis on statistical errors will no longer be the sole driving force. The understanding, control, and reduction of systematic errors is the next major challenge—there are several cases where measurements are already systematics-limited. Additionally, in the nonlinear regime of structure formation, state of the art measurements are pushing up against the limits of available prediction accuracy from theoretical computations and associated modeling. Major simulation campaigns, in concert with observational inputs, will be needed to control the error budgets. The simulation-based theoretical modeling process—consisting of a combination of first-principles calculations and parameterized models—will have to account accurately, and robustly, for the nonlinear evolution of structure formation, gas physics, and complex feedback effects.

Two examples will suffice to provide illustrations of the sort of difficulties that have to be faced. The first is cluster cosmology—it was recently shown by Cunha & Evrard (2010) and Wu et al. (2010) that the halo mass function has to be predicted at the 1% level of accuracy for a DES-like survey in order to avoid errors on dark energy parameters. At the same time, an error in the prediction of halo masses at the 2% level will translate into mass function uncertainties in the cluster mass range of 5%–10% (see Bhattacharya et al. 2011). The effect is particularly strong because of the exponential fall-off of the mass function at high masses. The second example relates to prediction error requirements for weak lensing shear. These are also severe—for the matter power spectrum, predictions at the percent level accuracy at k ∼ 10h⁻¹ Mpc are needed to extract all available information on dark energy (Huterer & Takada 2005). A more recent study by Hearin et al. (2012) concludes that the accuracy requirements are even more stringent—for imaging surveys such as LSST and Euclid, the power spectrum needs to be predicted at 0.5% accuracy out to k ∼ 5 h⁻¹ Mpc.

Currently, it is common practice to use fitting functions motivated by theoretical heuristics and matched to simulations to provide the required modeling input. For instance, analyses of weak lensing measurements employ fitting functions for the matter power spectrum. Semboloni et al. (2006) used the Peacock and Dodds fitting function (Peacock & Dodds 1996) as well as Halofit (Smith et al. 2003) for their analysis of the CFHTLS cosmic shear measurements. Lin et al. (2012) used Halofit for their recent SDSS analysis, as do Kilbinger et al. (2013) in the most recent CFHTLens analysis (they point out that an emulator would be more accurate but one was not available over the required redshift range) and Jee et al. (2013) for the DLS analysis. Over a limited k range, the fitting functions are accurate at the ∼5%–10% level for ΛCDM models (see, e.g., Heitmann et al. 2010). The latest improvement of Halofit is provided by Takahashi et al. 2012; a comparison of our results with theirs is provided in Section 6.3.

For next-generation applications, the fitting function approach suffers from two deficiencies, first, error-control is non-uniform over cosmological model parameter space, and, second, the fitting process degrades the actual accuracy of the simulations used to build the parameterized fit. Moreover, because the fitting forms become essentially arbitrary as the accuracy requirements become more stringent, this strategy is difficult—if not impossible—to implement systematically across a large number of freely floating cosmological and modeling parameters.

Cluster mass functions provide another example of these difficulties—obtaining dark energy constraints from the abundance of clusters of galaxies as carried out in, e.g., Vikhlinin et al. (2009), employs fitting forms for the mass function; to extend these fits across cosmological models, the assumption of universality (Jenkins et al. 2001) is required. However, as shown in, e.g., Bhattacharya et al. (2011), universality for wCDM models only holds at the ∼10% level of accuracy. It appears quite unlikely that a simple mass function fit—valid for a large class of dark energy models, covering a wide redshift range, and satisfying percent level accuracy requirements—will ever be attainable.

We have recently developed the "Cosmic Calibration Framework" (CCF) to provide accurate prediction schemes for cosmological observables (Heitmann et al. 2006, Habib et al. 2007) seeking to avoid the shortcomings of the fitting function approach mentioned above. The CCF contains error-controlled direct numerical oracles for the predicted quantities, as opposed to (potentially uncontrolled) analytic fitting forms. It aims to make tools available that provide essentially instantaneous predictions of large scale structure observables, e.g., the nonlinear power spectrum, mass functions for different halo definitions, or the halo concentration–mass relation.

At the heart of the CCF lies a sophisticated sampling scheme for optimally placing a given finite number of cosmological models in parameter space. Simulations are carried out at these parameter values (we use orthogonal-array Latin hybercubes as well as symmetric Latin hybercube designs). The next step is an efficient representation that translates the measurements from the simulations into functions that are conveniently interpolated (we use a principal component basis to provide a reduced data representation), and finally an accurate interpolation scheme over the basis functions (our choice here is Gaussian process modeling). For an introduction to the general framework, see, e.g., Santner et al. 2003.

The CCF was first described in Heitmann et al. (2006), with details and examples given in Habib et al. (2007). In a series of three papers (Coyote Universe I–III), we described the development of a precision emulator for the matter power spectrum over a five-dimensional parameter space (θ = {ω_b, ω_m, n_s, w, σ₈}). This emulator provides predictions for the power spectrum for wCDM cosmologies out to k ∼ 1 Mpc⁻¹ at the 1% accuracy level for a redshift range of 0 ⩽ z ⩽ 1, over a relatively wide range of cosmological parameters (see Section 3). Since then, similar approaches have been followed in Schneider et al. (2011) and Agarwal et al. (2012) with some differences in the sampling and interpolation schemes used to build power spectrum emulators.

The CCF framework was extended in Schneider et al. (2008) to derive an approximate statistical model for the sample variance distribution of the nonlinear matter power spectrum. Eifler (2011) used the emulator to generate a weak-lensing prediction code to calculate various second-order cosmic shear statistics, e.g., shear power spectrum, shear-shear correlation function, ring statistics and Complete Orthogonal Set of EB-mode Integrals (COSEBIs). The emulator we presented in Lawrence et al. (2010) was used by Huff et al. (2011) in their analysis of SDSS weak lensing (where it was combined with Halofit to obtain the power spectrum at larger k values). More CCF-based emulators are under development; a recent example is an emulator for the halo concentration–mass relation (Kwan et al. 2013).

Observational requirements, such as those for DES weak lensing, and the work by Eifler (2011) and Huff et al. (2011) have prompted us to extend the Lawrence et al. (2010) power spectrum emulator in three directions: (1) making the Hubble parameter h freely choosable (in the original emulator the value for h in each model is fixed by the CMB distance to last scattering constraint); (2) extending the k range out to k = 8.6 Mpc⁻¹; (3) extending the z range out to z = 4. As described in Section 3 below, these extensions are nontrivial—a brute force approach would require simulations of size L = 1.5–2 Gpc with at least $N_p^3=10,000^3$ particles to fulfill the requirements for 1% accuracy, as derived in Heitmann et al. (2010). The most stringent requirement is due to the fact that the power spectrum should not be measured further out than half of the particle Nyquist frequency, k_Ny = π/Δ_p with Δ_p being the initial particle separation L/N_p, and at the same time the linear modes should be well resolved in the simulation volume. In addition, at high redshift and large k, particle shot noise becomes a significant issue (see the discussion in Section 3).

Running a large number of the brute force simulations described above is currently impractical. In order to extend our prediction scheme we therefore choose another route, by using a set of nested simulations of various volumes, and matching the results together at different values of k and z. The box lengths are chosen to be L = 1300 Mpc (these are the simulations from the original Coyote Universe project as described in Lawrence et al. 2010) and L = 365 Mpc below z ∼ 0.7. For higher redshifts and smaller scales we use simulations of size L = 180 Mpc, and for redshifts z > 2, L = 90 Mpc.⁶ The specific choices for the box size are explained in more detail in Section 3. In short, these choices ensure good overlap between the resulting power spectra at different scales and redshifts and limit inaccuracies due to shot noise and finite box size effects.

A nested approach in the same spirit, and with similar matching scales in k and redshift, has been successfully employed by Takahashi et al. (2012). Based on their simulation results, these authors have developed an improved version of Halofit for predicting the nonlinear power spectrum. We compare their results to ours in detail in Section 6.3. The nested box approach has its deficiencies, as discussed further in Section 3. Nevertheless, carrying out a battery of tests we can demonstrate that our predictions hold at the level of better than 5% error over the full k range. Appendix A shows one example of a large, high resolution simulation that spans most of the k-range covered by the new emulator and the comparative error is ∼±2%.

At small length scales, complex baryonic processes become important. Because these are difficult to model, a relatively large uncertainty in the power spectrum exists over the upper k range treated in this paper. (Massive neutrinos also have a small effect, see, e.g., Bird et al. 2012.) White (2004) estimated the effect of baryonic cooling on the mass power spectrum using the halo model, finding an increase in power at k ∼ 10 h Mpc⁻¹ at the level of a few percent at z = 0. An analogous approach by Zhan & Knox (2004) to study the effect of the intracluster medium (hot baryons) found a suppression in the power spectrum at roughly similar levels.

Shortly after these papers appeared, simulations were carried out to study these effects in more detail. Jing et al. (2006) carried out two simulations with the smoothed-particle hydrodynamics (SPH) code Gadget-2, a non-radiative gas simulation and another including radiative cooling and star formation, with a third gravity-only simulation serving as a reference. At k ≃ 10 h Mpc⁻¹ they found an effect on the total matter power spectrum at the few percent level in the non-radiative gas simulation and at the 10% level in the gas simulation with radiative cooling and star formation at z = 0 (the effects are smaller at higher redshifts—for the radiative cooling/star formation simulation, they are at the 2%–3% level at z = 1). In both cases, the baryonic effects led to an enhancement in the power spectrum. A similar campaign was later carried out by Rudd et al. (2008): gravity-only, non-radiative gas dynamics, and a third simulation including radiative heating and cooling of baryons, star formation, and supernova feedback. This set of simulations used the Adaptive Refinement Tree (ART) code, combining an N-body method and an Eulerian hydrodynamics solver on an adaptive mesh. For the non-radiative gas simulation they found a similar effect as Jing et al. (2006) on the overall matter power spectrum, though effects on the baryon and dark matter components separately were different. The results from their third simulation showed a dramatic increase of the matter power spectrum well beyond the effect that Jing et al. (2006) observed (at the 50% level at k ≃ 5  h Mpc⁻¹)—but as noted by the authors, this result was not meant to be definitive.

More recent studies have been carried out by van Daalen et al. (2011), again using Gadget-2 (the authors also provide a much more comprehensive discussion of results from other groups than we have space to do here). The main difference in this study is the inclusion of feedback from active galactic nuclei. They observe a strong effect on the power spectrum in the opposite direction than found by the previous studies: a decrease in the matter power spectrum at 1% at k ≃ 0.3 h Mpc⁻¹, at 10% at k ≃ 1 h Mpc⁻¹, and at 30% at k ≃ 30 h Mpc⁻¹. As shown by these results, obtaining accurate first principles predictions for the power spectrum including baryonic effects will be very difficult. Multiple issues, ranging from physical effects to simulation uncertainties, remain to be sorted out.

A strategy to properly deal with these difficulties is still emerging, but it will definitely require methods to incorporate information from observations, baryonic simulations, and an accurate calibration of the gravity-only power spectrum (which has no free modeling parameters). The difficulty of accounting for all "gastrophysics" effects at the desired accuracies, almost certainly precludes using hydro simulations in a fully predictive mode. One possible pathway is to construct parameterized models for baryonic effects on top of predictive N-body simulation results (see e.g., Semboloni et al. 2011; Zentner et al. 2013 for some recent work) and combine these with (possibly multiple) observations ("self-calibration"), to successfully extract cosmological information. In this paper, our target is to establish the bedrock on which all such analyses must rest—an accurate calibration of the gravity-only matter power spectrum.

The rest of the paper is organized as follows. After providing a brief summary of our power spectrum estimation from the simulations (Section 2), we describe the cosmological model space covered and the simulation suite used to build the emulator (Section 3). In Section 4 we outline our strategy to include the Hubble parameter h as a free parameter, without adding more N-body simulations. Next we discuss the generation of smoothed predictions for the power spectra for each model; this process underlies the interpolation scheme for constructing the emulator. In Section 6 we show some examples from the working emulator, including a brief comparison with Halofit. Finally, we end with a conclusion and outlook in Section 7. Appendix B presents a short discussion of an improvement to the general Gaussian process approach employed here compared to that in the previous Coyote papers, including a minor numerical correction of the earlier results.

The extension of the original emulator beyond k = 1 Mpc⁻¹ and z = 1 at high accuracy is nontrivial due to two limiting factors: the spatial Nyquist frequency,

$$\begin{equation} k_{\rm Ny}=\frac{\pi N_p}{L}, \end{equation} \tag{ 6 }$$

setting the largest viable k-value, and the particle shot noise limit:

$$\begin{equation} P_{\rm shot}(k)=\left(\frac{L}{N_p}\right)^3, \end{equation} \tag{ 7 }$$

restricting the lowest amplitudes at which the power spectrum can be measured accurately. These issues have been known for a long time (see, e.g., Baugh et al. 1995 for an early discussion). With the recent efforts to obtain very accurate results for the absolute measurement of the power spectrum out to small scales, it is worthwhile to briefly revisit the essential arguments.

Both the Nyquist limit and the shot noise level depend on the mass resolution and the size of the simulation volume and in both cases reasonable results can only be obtained for a large number of particles. Computational limits imply a necessarily finite number of particles; the obvious option is to shrink the simulation volume. But this option has its own pitfalls—lack of long-range power and increased sampling variance, for example—and will break down at some point as discussed quantitatively below.

The shot noise effect is particularly annoying at early times; at late times, once the power spectrum has risen substantially above the shot noise level, it becomes much less problematic. The Nyquist sampling issue, which essentially sets the dynamic range at which the initial condition can be produced, is worse for larger boxes with the same particle number. In addition, the impact of both problems depends on σ₈. For two simulations that differ only in their value of σ₈, the one with the higher σ₈ will contain more nonlinear structure and growth at the same redshift—this means that the shot noise problem is worse for low σ₈ runs while the issues with small volumes at lower redshift are more severe for high σ₈ models.

At the current levels of available computational power, it is impossible to carry out a brute force approach to obtain a sub-percent accurate answer for the power spectrum at the desired k-values. As discussed at length in Heitmann et al. (2010), the maximum value is set by k_Ny/2. At the same time, in order to ensure that the largest modes in the box at z = 0 are linear to high accuracy, a linear box size of 2000 Mpc would be a good choice. To reach a wavenumber value of 10 Mpc⁻¹ would imply a particle loading of $N_p^3\sim 12700^3\sim 2$ trillion particles. While it is possible to carry out a few such simulations (Habib et al. 2012), a large suite of them is clearly currently out of reach.

This notional set-up would lead to a shot noise level of ∼4 × 10⁻⁴. The model with the lowest σ₈ in our simulation set is M019 with σ₈ = 0.6159. The amplitude of the power spectrum at z = 4 is P(k = 10, z = 4) = 0.025 Mpc³, translating to ∼60 times above the shot noise level and is, hence, safe. If we wanted to decrease the number of particles by choosing a smaller volume, 1000 Mpc would be barely large enough (for more detailed discussions of finite volume effects see e.g., Heitmann et al. 2010). In this case, ∼250 billion particles would be needed to push out the Nyquist limit far enough. This would lead to a shot noise level of ∼4 × 10⁻³, dangerously only a factor of 6.25 below the amplitude of the power spectrum at k = 10 Mpc⁻¹ and z = 4. Moreover, even this simulation size constitutes a currently prohibitive expense if a full suite of simulations has to be performed.

Due to these obstacles, a strategy for mixing and matching simulation sizes—differently for different redshifts—needs to be employed. Smaller simulation volumes are required for high redshift results while larger volumes will be used for lower redshifts. Although such a procedure is effective as we show below, because of the uncertainties induced by having to sew together multiple simulation results, and because the number of cosmological models used is still the same as in the original set, sub-percent accurate power spectra are not obtained over the new, much wider dynamic range in wavenumber and redshift. At the small spatial scales that the current emulator goes to, however, current uncertainties in characterizing baryonic effects clearly outweigh its inaccuracies. Therefore, over the extended dynamic range, the constraints on the accuracy of the dark matter power spectrum are not uniform, and not as stringent near the limits of the range probed in wavenumber and in redshift.

We quantify below some of the errors due to the fact that the power spectra at all k and z do not arise from a single overarching simulation. Over its full range, and across all cosmological models, the emulator is nevertheless accurate to better than 5%. We now describe the details of our matching strategy for the power spectra at each redshift. We store results at 11 scale factors (a = 0.2, 0.25, 0.2857, 0.3333, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, see also Tables 3–5) and then interpolate between those outputs to obtain results at any scale factor in between.

For the largest scales, we use renormalized perturbation theory (RPT; see Crocce & Scoccimarro 2006 and Crocce & Scoccimarro 2008 for details on the underlying idea) using the Copter code (Carlson et al. 2009); the code has been modified to allow for w ≠ −1. We start by breaking up the models into three groups, depending on the power spectrum normalization, σ₈:

$$\begin{eqnarray} &&\sigma _8 < 0.7, {\rm Assembly\; 1}, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &&0.7\le \sigma _8\le 0.8, {\rm Assembly\; 2}, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &&\sigma _8 > 0.8, {\rm Assembly\; 3}. \end{eqnarray} \tag{ 10 }$$

We use RPT for all three cases in the same way up to the following k values:

$$\begin{eqnarray*} &&1.0 \le a\le 0.5: k=0.03 \,{\rm Mpc}^{-1}\nonumber \\ &&0.4 \le a\le 0.2: k=0.1 \,{\rm Mpc}^{-1}.\nonumber \end{eqnarray*}$$

(As described above, we do not store any outputs between a = 0.4 and a = 0.5.) The first cut is the same as was used in Lawrence et al. (2010), while the second is more aggressive. Since the power spectrum is more linear at higher redshifts, perturbation theory will be valid out to higher k. We exhibit the accuracy of RPT at the matching scales for the model with the highest σ₈, M037 (σ₈ = 0.9)—the most difficult case—in Figure 1. The ratio of RPT is shown with respect to the simulation output (the average from our 20 realizations with L = 1300 Mpc) at three different redshifts, z = 0, 1.5, 4. The red vertical line shows our matching point for perturbation theory. In all cases, the error is below 1%, in very good agreement with the findings of Carlson et al. (2009).

Zoom In Zoom Out Reset image size

Next, we must determine the maximum k values out to which to use the results from each of the different volumes. The individual matching strategies for the three assemblies depending on σ₈ are summarized in Tables 3–5. As mentioned above, the matching point depends on the Nyquist wavenumber and, at higher redshifts, the shot noise level. As previously established in Heitmann et al. (2010) the results from the 1300 Mpc volume hold at the sub-percent level accuracy out to k ∼ 1 Mpc⁻¹ between 0 ⩽ z ⩽ 1 or 0.5 ⩽ a ⩽ 1.0. Therefore, the power spectra for all three assemblies in this range are taken from these simulations (the simulations are also particularly valuable because of the large numbers of realizations for each of the boxes).

At higher redshifts, the shot noise level of the earlier simulations becomes too high and we can use them only to a lower k-value, consequently, in all three assemblies, for 0.2 ⩽ a ⩽ 0.4 we cut off the power spectra from the 1300 Mpc boxes at k = 0.4 Mpc⁻¹. For a ⩾ 0.5, the upper k value is taken to be 1.0 Mpc⁻¹. The following step is to match the results from the smaller boxes at these cutoff points. Shot noise level considerations play the dominant role in choosing the matching points for the power spectra from the different simulations.

In the case of Assembly 1 (high σ₈), the second-largest box (365 Mpc) is used to cover the remaining full range in wavenumber from 0.4 ⩽ a ⩽ 1. The next two smaller boxes of size 180 Mpc and 90 Mpc are used to fill in the higher k range as shown in Table 3. The results of analogous strategies for Assembly 2 (medium σ₈) and Assembly 3 (low σ₈) are given in Table 3 and Table 5.

Figure 2 shows the matched power spectra for two example models, M019 and M037, representing the lowest and highest σ₈ in our model selection. The upper panels show the full power spectra from all boxes while the lower panels show the results with cut-offs applied at the appropriate k-values. The results in the lower panel demonstrate that the approach of using nested volumes is qualitatively satisfactory. In the next two sub-sections we discuss the main errors involved in the matching procedure: (1) finite mass resolution which leads us to push beyond k_Ny/2 for the small boxes, and (2) finite volume effects.

Zoom In Zoom Out Reset image size

In Heitmann et al. (2010) the effects of finite mass resolution on power spectrum estimation were discussed in some detail. Two effects investigated there act in opposite directions. The first is sampling noise, i.e., the effect on computing a power spectrum from a density field computed from a different number of particles. For this, a particle distribution at z = 0 was down-sampled by factors of 8, 64, and 512 and the power spectra re-measured. The effect of the insufficient sampling of the density field led to an enhancement of the power spectrum larger than 1% beyond k_Ny/2. The second is the effect of lower particle sampling in the initial conditions. This leads to a suppression of the power spectrum since the initial conditions now lack small-scale power. This effect decreases for simulations with larger Nyquist frequency since the number of modes in the box at small scales increases. In the test in Heitmann et al. (2010) it was shown that a simulation with k_Ny = 0.859 h Mpc⁻¹ has lost more than 10% of power at k_Ny/2, while for k_Ny = 3.44 h Mpc⁻¹ it is well below 1% at k_Ny/2.

In the current paper, we use the large volume simulations (1300 Mpc) only up to k ⩽ 1 Mpc⁻¹ at low redshift and k ⩽ 0.4 Mpc⁻¹ at high redshift, both values being well below k_Ny/2. For the smallest volumes (90 Mpc) we also do not use any results beyond k_Ny/2, simply because for the smallest boxes the Nyquist limit is much larger than k ∼ 10 Mpc⁻¹. For the intermediate boxes (180 Mpc and 365 Mpc) we choose to be more aggressive. For the lowest redshifts we push the 365 Mpc limit to twice the Nyquist criterion and for intermediate redshifts we use the 180 Mpc box out to the Nyquist frequency. In principle, we could have avoided pushing the 365 Mpc box beyond the Nyquist limit by using the 180 Mpc box at large k values instead. There are two reasons why we chose not to do this.

First, the comparison of inaccuracies due to finite box effects versus going beyond the Nyquist limit suggests that the finite box errors dominate. The quasi-linear regime in the small box simulations close to the final time step at which they are used is not accurate because of the missing large-scale power. In addition, the scatter in the overall amplitude is larger due to the small box size as discussed in the next sub-section. These two errors combine to outweigh the inaccuracy due to the Nyquist limit. This point is illustrated in Figure 2. In the right upper panel (M037), at z = 0, the green curve shows the result for the 365 Mpc box while the blue curve shows the result for the 180 Mpc box. It is apparent that the small box clearly does not capture the nonlinear turn-over behavior beyond k ∼ 0.1 Mpc⁻¹ very well. On the other hand, both results are in close agreement at large k, so the accuracy of the 365 Mpc box results is still good at that point.

Second, the matching procedure described in Section 5 for melding the different power spectra pieces induces another small inaccuracy. This is again mainly due to finite box size effects, as the amplitude has to be adjusted to enable a matching of the boxes in the high k region. It is therefore desirable to minimize the number of matching points, and hence to take the results from the 365 Mpc boxes out to higher k values, rather than to introduce another matching point at intermediate scales. In order to verify that using the 365 Mpc box at higher k does not introduce a major error, we also checked that the difference at high k between the 365 Mpc box and the 180 Mpc box is small, as can be seen in Figure 2. For M037 (high σ₈) we find differences of 2% and less beyond k = 6 Mpc⁻¹ and for M019 (low σ₈) we find differences well below 1% in the high k range.

To summarize this discussion, for the high redshift results, we do not use simulations beyond k_Ny/2 (see Tables 3–5 for more details) and therefore estimate the error throughout the k range as being well below 1%, based on the test results from Heitmann et al. (2010). For the low redshift cases, we use some results that go beyond the Nyquist limit. By comparing to the smaller boxes we estimate that the error does not exceed 2%.

As outlined in the Introduction, ultimately one would want to carry out simulations in large box sizes, 1–2 Gpc, but this is currently impractical for a large suite of runs, each with sufficient mass resolution. Smaller simulation volumes have two drawbacks: (1) At redshifts close to z = 0 the undersampling of large-scale power may be a problem (the precise extent depending on σ₈ and the box size). In smaller boxes, while the largest-scale modes might appear to be linear, their amplitude is actually suppressed by missing nonlinear power (in a large volume run, the same modes would have higher amplitudes). (2) The realization scatter in small volumes is much larger due to the smaller number of large scale modes. The first of these points was considered in depth in Heitmann et al. (2010). There, we demonstrated the suppression of the power spectrum due to finite volume effects by comparing the average of 4 realizations in a 2000 h⁻¹ Mpc box, 8 realizations in a 936 h⁻¹ Mpc box, and 127 realizations in a 234 h⁻¹ Mpc box. Figure 6 in Heitmann et al. (2010) shows the suppression at mildly nonlinear scales. The realization scatter was not a significant issue in that paper since the final simulation volumes were large (1300 Mpc) and we averaged over 20 realizations.

In the current paper, the realization scatter is of concern since we do resort to using small volumes. The mildly nonlinear regime is still very accurate since it is covered by large volume simulations, so the key question is to understand the high-k behavior. In order to investigate this effect, we carry out two tests. First, we use the same 127 simulations used in Heitmann et al. (2010) to measure the overall dispersion in the power spectra between different realizations. These simulations were carried out in 325 Mpc volumes with a PM code almost matching the small simulation volumes used here. Since we are only interested in the relative effect, these lower resolution simulations are sufficient. The results are encapsulated in Figure 3. The upper panel shows the ratio of the average of the 127 power spectra with respect to each individual spectrum; in the lower panel, we attempt to correct the run-to-run scatter by simply adjusting the amplitude of each power spectrum to the average value at the highest measured k (which is well-sampled in each simulation). This procedure improves the scatter somewhat but is rather uncontrolled and we decided not to use it. Nevertheless, the result is somewhat interesting. Overall, the run-to-run scatter at small scales is up to about 5%, at the matching scale of k ∼ 1 Mpc⁻¹, and up to 10% for the most extreme outliers. In order to reduce this problem, we carry out at least two realizations for each model, in some case up to four to obtain a realization in which the matching to the 1300 Mpc box is as close as possible. Since we match the even smaller volumes (180 Mpc and 90 Mpc) at larger k values, the problem there is not quite as severe. Nevertheless, even in those cases we run up to four realizations per model to provide an accurate match to the larger boxes.

Zoom In Zoom Out Reset image size

As a second check we carry out six simulations of the same cosmology each with box length L = 365 Mpc and $N_p^3 = 512^3$ particles, which are chosen to exactly match the size and resolution of the smaller Coyote runs. The initial particle distribution is set with the same realization at z = 211 on large scales for all six simulations, but on scales smaller than k = 1.39 Mpc⁻¹, we allow the initial conditions to vary between runs by changing the random seed for each simulation. These are then evolved using Gadget-2 and we show the ratio of power spectra at z = 0 with reference to a single simulation in the set in Figure 4. The red line marks where the initial power spectrum differs between the realizations. Notice that the scatter in the power spectrum is less than ∼3% and appears unbiased. This gives an estimate of the amount of mismatch between the 1300 Mpc and 365 Mpc Coyote simulations that we can expect after matching their power spectra. There is also a small amount of leakage of power from small to large scale modes, since the power spectra only match at k ≪ 1.39 Mpc⁻¹, where they were seeded identically in the initial conditions.

Zoom In Zoom Out Reset image size

In terms of the effects of the missing low-k power, at z = 0, for a ∼100 Mpc box, the rms amplitude of the DC mode is at the level of ∼10%, falling to a few percent at ∼300 Mpc (see, e.g., Gnedin et al. 2011). These numbers are sufficiently small—given the level of accuracy we are aiming to attain here—that a more sophisticated correction procedure is not required. As previously mentioned, box sizes in the 1–2 Gpc range are sufficiently large to render these effects sub-dominant when targeting accuracy levels of ∼1% Heitmann et al. (2010). As more computer power becomes available, construction of the next generation of emulators will profit from the increased volume and better mass resolution, both of which are essential to improve the accuracy at higher wavenumbers.

To summarize, finite box size effects are clearly an important issue. Some of the problems such as realization scatter can be overcome by generating a large number of realizations (though this is expensive) but some small inaccuracies in the power spectrum will be unavoidable until larger volume, high-particle loading simulations can be carried out.

As explained in Section 2 the emulator is built using results from 37 cosmological models, specified in Table 1. In addition to these models, we have generated P(k) results for a ΛCDM model (M000), which is close to the current best-fit measurements. Simulations for this model were carried out in exactly the same way as for the other models, using several realizations for each box size, and then matching them to produce a smoothed power spectrum. This reference P(k) can be used as one test of the emulator's accuracy.

Figure 8 shows results for two emulator versions, one in which the value of h is locked by the CMB distance to last scattering constraint for each model, as in the original emulator (Lawrence et al. 2010) (upper panel), and a second version where h is allowed to be a free parameter following the approach presented in Section 4 (lower panel). Both panels show the ratio of the emulator to the simulation result as a function of k, the different solid curves correspond to eleven different redshifts between z = 4 and z = 0 used to build the emulator. In the h-fixed case, over a wide k range, and for all redshifts, the emulator prediction is accurate at the 1% level, degrading slightly only at larger k values, but still remaining better than 5%.

Zoom In Zoom Out Reset image size

When h is allowed to range freely (lower panel), then over the k range where RPT is used (up to k ∼ 0.03 Mpc⁻¹ for the low redshift results), the result is excellent, better than that in the top panel. This is not surprising since for this range we now have 100 models to predict the power spectrum. Beyond that matching point, the accuracy degrades slightly in the quasi-linear regime compared to the more restricted emulator, but is still around 1%. Finally, beyond k ∼ 1 Mpc⁻¹, the predictions are less accurate, but the worst case error remains around 5%.

Next, we present results from a number of "holdout" tests. In these tests, one model out of the 37 is excluded and an emulator is constructed based on the remaining 36 models. Then with this new emulator, a prediction for the excluded model is generated and the accuracy of the prediction determined. This test has an obvious caveat, particularly relevant in cases where not many models are available—degradation of the emulator quality because an important point in the parameter-space hypercube is omitted. This can be particularly serious if the excluded point is on the edge of the design, because now extrapolation to the edge of the hypercube is required, which has a much higher error. In order to avoid this additional complication, we restrict our holdout test to models that lie well within the hypercube. Nevertheless, it should be kept in mind that reported errors in the holdout tests can be larger than the actual values. The holdout test results, shown in Figure 9 are similar to the ones found for M000, both in error amplitude and trends with increasing k. In the case where h is allowed to be a free variable, the predictions are again less accurate, already in the quasi-linear regime, and degrade further beyond k ∼ 1 Mpc⁻¹, to of order 5%. Note that the low-k error band is consistent with that for the linear theory test case (see Figure 5).

Zoom In Zoom Out Reset image size

To summarize, the accuracy of the new, extended, emulator is well-characterized at all redshifts. With h fixed to the CMB constraint, the accuracy to k ∼ 1 Mpc⁻¹ is at the ∼1% level, degrading to ∼5% for k ∼ 10 Mpc⁻¹. When h is allowed to be a free parameter, the degradation in accuracy is relatively modest.

In order to improve the accuracy of the emulator with h included as a free parameter, a new design with more data points would have to be created. Based on the convergence tests in Habib et al. (2007), where the design space was varied between 32 and 128 models, a design with around 40 to 50 cosmological models should lead to percent level errors. We did not follow this strategy in the current paper because it would not allow us to use the simulation data already obtained in our previous work. Some of these questions (e.g., how to create nested designs that allow improvements of the emulator accuracy by adding more models in a systematic fashion) are currently under investigation, and will be addressed in future work.

With the full emulator at hand we can now carry out a sensitivity analysis for the nonlinear power spectrum, along the same lines as conducted in Section 4 for the linear power spectrum. As before, we fix five of the six cosmological parameters at the midpoints of their range and then vary the sixth parameter between its minimum and maximum value. The results—for the lowest and highest redshift we consider—are shown in Figure 10. The upper six panels show the results at z = 4 and the lower six panels at z = 0.

Zoom In Zoom Out Reset image size

Some interesting features in the results can be noted. As to be expected, the baryon fraction ω_b does not affect the power spectrum significantly. The best constraints we have for ω_b come from CMB measurements; Planck determines this parameter to exquisite accuracy. The effect of the matter fraction ω_m dominates at large scales, here again, CMB results deliver very good estimates of this parameter. In general, it is important to allow for variations in ω_m due to degeneracy issues; ω_m influences the overall amplitude of the power spectrum as do σ₈ and h. The spectral index n_s, influencing the tilt of the power spectrum also mainly effects large scale behavior. The remaining three parameters, h, w, and in particular σ₈, alter the power spectrum on all scales and constraints on these parameters can therefore be improved by measurements of nonlinear scales.

In particular, the influence of σ₈ at the two redshifts shown in Figure 10 on the power spectrum is noteworthy. At high redshift, nonlinearities for models with large σ₈ (pink) have already developed considerably, therefore the ratio with the midpoint power spectrum shows large values. The power spectra for small values of σ₈ (blue) show only mild nonlinear growth and therefore the ratio is still almost flat. Over time, the nonlinear turn-over moves in to smaller k ranges (affecting larger and larger scales), causing the bump just before k = 1 Mpc⁻¹ shown in the panel for z = 0. At this epoch the overall difference between the models has decreased, since at this point all models have reached the nonlinear regime. A similar observation can be made for w—at early times the differences between the models are more pronounced than at later times. It is easy to imagine using the new emulator as a convenient tool to investigate the effects of varying cosmological parameters and redshift on P(k) in the nonlinear regime.

We now perform a comparison between the emulator and Halofit, a popular fitting formula for the matter power spectrum motivated by the halo model. We use the updated version with 35 fitting parameters provided by Takahashi et al. (2012), which is calibrated to larger volume and higher resolution N-body simulations than the original Smith et al. (2003) formula. It is based on simulations run on six WMAP cosmologies (years 1, 3, 5, 7, and two WMAP7 models, but with w = −0.8 and w = −1.2, instead of w = −1). The obtained fitting formula was checked against simulations run using the first 10 of the 37 Coyote models. This improved version of Halofit is stated to be accurate to ∼5% in the range 0 < k < 1 h Mpc⁻¹ and ∼10% at 1 < k < 10 h Mpc⁻¹, for 0 < z < 3.

One important feature of the emulator-based approach is that, close to the parameter values of the underlying simulations, the emulator error is small, if not essentially the simulation error. However, this is not the case with a fitting formula. For instance, even though the N-body simulations that underlie the new Halofit formula agree with the original Coyote emulator at the ∼1%–2% level over its range of validity for the 10 Coyote models used as a reference in Takahashi et al. (2012, see Figure 2), this good level of agreement drops to 8%–13% (as a function of the k range) when the fitting formula is used in the comparison. Depending on what one may be trying to achieve, this loss of accuracy in Halofit is significant.

We first compare the predictions of the emulator (both versions, h free and h fixed) at z = 0 for M000 against the smooth simulations (in the same way as shown in Figure 8) and the prediction of Halofit for the same model. Since the emulator was built without including this simulation, this is a good test of both prediction schemes, emulator and Halofit. The results are shown in Figure 11. The emulator predictions for the h-fixed version are accurate at 1% out to k ∼ 1 Mpc⁻¹ and then degrade very slightly to 2%. The h-free version of the emulator is accurate at the 1%–2% level throughout. Halofit shows 3% deviations compared to the simulation in the mildly nonlinear regime (k ∼ 0.1 Mpc⁻¹) and more than 6% in the nonlinear regime.

Zoom In Zoom Out Reset image size

Next we show a comparison of the emulator (h-free version only) directly with Halofit for a variety of cosmologies at z = 0. Figure 12 displays the ratio of the emulator with respect to linear theory power spectra from CAMB and Halofit as modified by Takahashi et al. (2012). We choose to test a WMAP7 (Komatsu et al. 2011) and a Planck (Ade et al. 2013) cosmology and two other cosmologies on the limits of the design of our emulator. The red and blue curves in particular are on the edge of the new design with h as a free parameter. The cosmological parameters for these two model are the same as model M000: ω_m = 0.1296, ω_b = 0.0224, n_s = 0.97, σ₈ = 0.8, w = −1, h = 0.72 and in our comparison we vary one to two parameters at a time as noted on the legend.

Zoom In Zoom Out Reset image size

At large scales, k < 0.01 h Mpc⁻¹, all power spectra are consistent at the percent level; this is less trivial than it seems since Halofit and linear theory start diverging at k ∼ 0.003 h Mpc⁻¹ and we use the full version of the emulator with extreme values of h located at the edge of the design. On quasi-linear scales, there is ∼5% deviation with some oscillatory structure from the BAO feature. On smaller scales, the accuracy is consistent with what is stated in Takahashi et al. (2012), apart from the matter dominated case with Ω_m = 0.51, in which the Halofit power spectrum is suppressed by more than 20%. The relatively uniform error behavior of an emulator—a result of the sampling theory and the Gaussian process based-regression methodology—is hard to reproduce in a fitting formula; the last result of the comparison provides an example of this.