THE COYOTE UNIVERSE. III. SIMULATION SUITE AND PRECISION EMULATOR FOR THE NONLINEAR MATTER POWER SPECTRUM

During the last three decades, cosmology has made tremendous progress: from order of magnitude estimates to measurements of key cosmological parameters approaching percent level accuracy. The standard model of cosmology, based on the growth of structure by gravitational instability, has been impressively validated on large scales where the perturbations to the Friedman model are small. While the model is successful at predicting or reproducing a wide array of observations, it contains several mysterious elements with perhaps the most mysterious being the accelerated expansion of the universe (Riess et al. 1998; Perlmutter et al. 1999). The accelerated expansion may be caused by a dark energy or hinting at a modification of general relativity on the largest scales.

To date, most of the best known cosmological parameters have been constrained primarily from the study of anisotropies in the cosmic microwave background (CMB) radiation. However, large-scale structure probes play an important role in breaking parameter degeneracies and in constraining conditions in the late-time universe. Such probes are becoming ever more precise, with future surveys aiming at measurements approaching percent level accuracy to better characterize the universe in which we live. Techniques based on the observation and analysis of cosmic structure include the use of baryon acoustic oscillations (BAOs), redshift space distortions, weak-lensing measurements, and the abundance of clusters of galaxies; they stand to play a pivotal role in improving our understanding of the dynamics of the universe. (For a recent discussion regarding improvements on dark energy constraints from combining different probes, see, e.g., Albrecht et al. 2006, 2009.) The large-scale structure of the universe contains information about both the geometry as well as the dynamics of structure formation. In combination, these two pieces of information can help distinguish between dark energy or a modification of general relativity as the prime cause of cosmic acceleration.

On small scales, the cosmological interpretation of structure formation probes is complicated due to the nonlinear physics involved. Commonly used fitting functions for, e.g., the power spectrum (Peacock & Dodds 1996; Smith et al. 2003) have poorly characterized systematics and are no longer adequate for precision work. Absent a controlled theoretical framework, direct use of simulations (augmented with phenomenological parameters as appropriate) is essential if the physics is to be more correctly captured. The simulation codes need to be adequately tested to ensure that they meet the new demands being placed upon them. The simulations which pass these goals are often very expensive and only a restricted number of runs can be performed. This in turn puts a premium on developing very efficient strategies to constrain parameters from limited observations and simulations. High-precision prediction schemes for different statistics are essential to succeed in this task.

This paper is the third in a series aimed at addressing this question in the context of the nonlinear matter power spectrum on Mpc scales. Current observations in weak lensing are quickly becoming theory limited due to the lack of precise theoretical estimates of this quantity for a wide range of cosmological models. This is a pressing problem. It is also relatively simple, allowing us to work through—in a concrete setting—the many steps that will be routinely required in the future. In some ways, the problem is one of the simplest currently confronting theorists, but as we shall discuss below even it has demanded collaboration with other communities, the development of a significant infrastructure, new modes of working, and large amounts of manpower and computational capacity.

In Paper I of this series (Heitmann et al. 2009b), we demonstrated that it is possible to obtain nonlinear matter power spectra with percent level accuracy out to k ≃ 1 h Mpc⁻¹ and derived a set of requirements for such simulations. Paper II (Heitmann et al. 2009a) described the construction of an emulation scheme to predict the nonlinear matter power spectrum and the underlying cosmological models, building on the "Cosmic Calibration Framework" (Heitmann et al. 2006; Habib et al. 2007; Schneider et al. 2008). Here, in Paper III, we present results from the complete simulation suite based on the cosmologies presented in the second paper and publicly release a precision power spectrum emulator.⁷ The simulation suite is called the "Coyote Universe" after the cluster it has been carried out on. We will extend our work to include other measurements, such as the mass function or higher order statistics, in future publications.

The outline of this paper is as follows. In Section 2, we describe the simulations we have run in support of this program and the codes with which they were run. Section 3 describes how we put together the different estimations of the power spectra from our multi-scale runs while Section 4 presents the details of the emulator and tests. We describe some lessons learned in Section 5 before concluding in Section 6.

2.1. The Simulations

2.2. The Cosmological Models

The selection of the cosmological model suite depends on two considerations: the statistical framework we use to construct the emulator and current parameter constraints from the CMB as set by WMAP-5 (Komatsu et al. 2009). We do not insist on a formal methodology to make the model selection, but instead apply some practical and conservative arguments to justify our decisions. For an in-depth discussion, we refer the reader to Paper II. In this paper, we briefly summarize some of the considerations since these will define the region for which the emulator will be valid.

Our aim is to find a distribution of the parameter settings—the simulation design—which provides optimal coverage of the parameter space, using only a limited number of sampling points. Simulation designs well suited for this task are Latin–Hypercube (LH)-based designs, a type of stratified sampling scheme. Latin hypercube sampling generalizes the Latin square for two variables, where only one sampling point can exist in each row and each column. A Latin hypercube sample—in arbitrary dimensions—consists of points stratified in each (axis-oriented) projection.

Very often LH designs are combined with other design strategies such as orthogonal array (OA)-based designs or are optimized in other ways, e.g., by symmetrizing them (more details below). By intelligently melding design strategies, different attributes of the individual sampling strategies can be combined to lead to improved designs, and shortcomings of specific designs can be eliminated. As a last step, optimization schemes are often applied to spread out the points evenly in a projected space. One such optimization scheme is based on minimizing the maximal distance between points in the parameter space, which will lead to more even coverage. Two design strategies well suited to cosmological applications in which the number of parameters is much less than the number of simulations that can be performed are optimal OA–LH design strategies and optimal symmetric LH design strategies. For this project, we have generated 40 different designs following different strategies and chosen the best design from these (where "best" refers to best coverage in parameter space with respect to specific distance criteria explained in Paper II).

In order to restrict the number of necessary simulation runs to be as small as possible, it is helpful to keep the number of cosmological parameters and their prior ranges small. Current observations of the CMB and the large-scale structure are consistent with a ΛCDM model with constant dark energy equation of state, w. We therefore concentrate on the following five cosmological parameters: ω_m ≡ Ω_mh², ω_b ≡ Ω_bh², n_s, w, and σ₈, where Ω_m contains the contributions from the dark matter and the baryons. We restrict ourselves to power-law models (no running of the spectral index) and to spatially flat models without massive neutrinos. For each cosmology, h is determined by the angular scale of the acoustic peaks in the CMB (Paper II) which is known to very high accuracy (0.3%; Komatsu et al. 2009). From WMAP five-year data, in combination with BAO, we have⁸

$$\begin{equation} \begin{array}{c} \omega _m = 0.1347\pm 0.0040\quad (3\%), \\ \omega _b = 0.0227\pm 0.0006\quad (3\%), \\ n_s = 0.9610\pm 0.0140\quad (2\%). \end{array} \end{equation} \tag{ 1 }$$

Current data restrict a constant equation of state for the dark energy to w = −1 with roughly 10% accuracy and recent determinations put the normalization in the range 0.7 < σ₈ < 0.9 with still rather large uncertainties.

Considering all these constraints and their uncertainties, we choose our sample space boundaries for the 37+1 models to lie within the range(s)

$$\begin{equation} \begin{array}{c} 0.120 < \omega _m < 0.155, \\ 0.0215 < \omega _b < 0.0235, \\ 0.85 < n_s < 1.05, \\ -1.30 < w <-0.70, \\ 0.61 < \sigma _8 < 0.9, \end{array} \end{equation} \tag{ 2 }$$

over which the emulator is designed to produce reliable results. We verified in Paper II that 37 models spanning these parameter ranges are indeed enough to generate an emulator at the 1% accuracy. We emphasize that our emulator is valid for the complete parameter space defined by these priors and not restricted to some values in a band around the best-fit cosmology (for current observational data). Obviously, the emulator quality will be slightly worse on the edges of the hypercube but the emulator is always accurate within ∼1% anywhere within the priors.

2.3. Power Spectra

In Paper I, we describe in detail how we obtain the matter power spectrum from a snapshot of the simulation. We briefly summarize the salient points here. We compute the dimensionless power spectrum,

$$\begin{equation} \Delta ^2(k) \equiv \frac{k^3P(k)}{2\pi ^2}, \end{equation} \tag{ 3 }$$

which is the contribution to the variance of the density perturbations per ln k. We obtain Δ² by binning the particles onto a 2048³ grid using CIC assignment, applying an FFT, correcting for the charge-assignment window function, and averaging the result in fine bins in |k| spaced linearly with width Δk ≃ 0.001 Mpc⁻¹. As discussed in Paper I, we do not correct for particle discreteness as our particle loading is high enough to make such corrections unnecessary and there are some indications that a simple Poisson shot-noise form is not correct (Paper I).

When creating a nonlinear power spectrum emulator, we prefer that the underlying training data set be smooth. We describe here how we construct a smooth power spectrum from the 20 realizations of each cosmology and from cosmological perturbation theory.

In each simulation, the modes in the initial density field are a single realization of a Gaussian random field and this introduces large run-to-run scatter. This scatter is reduced at higher k by the relatively large number of modes which are averaged. However, at low k the estimates of the power spectrum exhibit significant scatter which is expensive to reduce by brute force, i.e., running a very large number of realizations in large volumes. In addition, the approach to linear theory at low k can be quite slow for many currently popular models around ΛCDM (if 1% accuracy is the desired goal) so simply replacing the N-body results with the input linear model can be relatively inaccurate. This is however an area where perturbation theory can be of help, since the real-space mass power spectrum is computable in perturbation theory and there is a small but non-negligible range of scales where perturbation theory improves upon linear theory. For a recent overview of the performance of different perturbation theory approaches, see Carlson et al. (2009).

3.1. Perturbation Theory

While there has been a resurgence of interest in perturbative methods in recent years, we stick to the simplest and oldest method "standard perturbation theory" (Peebles 1980; Juszkiewicz 1981; Vishniac 1983; Goroff et al. 1986; Makino et al. 1992; Jain & Bertschinger 1994). We consider only the first correction to linear theory, which in standard perturbation theory can be written in terms of a simple integral over the linear theory power spectrum (we use the specific form given in Meiksin et al. 1999). When compared to our simulations, we find that standard perturbation theory is accurate at the percent level for k < 0.5 k_nl, where k_nl can be defined as (Matsubara 2008)

$$\begin{equation} k_{\rm nl}^{-2} = \frac{1}{3}\int \frac{dk}{k}\ \frac{\Delta ^2(k)}{k^2}. \end{equation} \tag{ 4 }$$

The values for k_nl at z = 0 and z = 1 are listed for all models in Table 1. For example, for model M000, k_nl ≃ 0.1 Mpc⁻¹ at z = 0. Similar results have been reported in Matsubara (2008) and Carlson et al. (2009).

Almost any scheme that switches smoothly from standard perturbation theory to our N-body results around 0.5 k_nl produces results that agree at the percent level. At z = 0, this would lead to a matching point between k = 0.045 Mpc⁻¹ (M036) and k = 0.075 Mpc⁻¹ (M019). For simplicity, we keep the matching point the same for all cosmologies. Since we have good statistics from our simulations at k ≃ 0.03 Mpc⁻¹ already, we choose this k value as a very conservative matching point for all models. To verify this point, we compare our simulation results to perturbation theory for the different cosmologies. The simplest approach for this comparison is "brute force": simply take the ratio of the simulation result to the perturbation theory prediction. The major obstacle here is the run-to-run scatter in the simulations. In order to overcome this problem, one can follow two routes: either incorporate fluctuations from the realization into the perturbation theory prediction, as done in Takahashi et al. (2008) or average over a large number of large volume simulations. We follow the second approach to avoid any ambiguities and systematic errors resulting from finite box size effects.

An example of our approach is shown in Figure 1 for a random subset of three models from the total of 37 used to build the emulator. Although at the lower end of the k range, the large but finite sampling volume becomes clearly evident, the matching to perturbation theory at k ≃ 0.03 Mpc⁻¹ works extremely well. A similar result can also be found in Carlson et al. (2009).

Zoom In Zoom Out Reset image size

3.2. The Estimation Procedure

Next, we show how we can combine perturbation theory and results from N-body simulations to generate smooth power spectra which will be the foundation for building our emulator. Two problems have to be solved for this. (1) We have to eliminate the scatter in the N-body results for the nonlinear power spectrum without erasing subtle features and match results from different resolution simulations. (2) We have to match very accurately between perturbation theory and simulation results. In the following, we discuss an approach based on process convolution to solve these problems.

3.2.1. Power Spectrum Estimation Using Process Convolution

In this section, we discuss the procedure for estimating the smooth power spectrum for each cosmology based on the simulation results which possess inherent scatter. Figure 2 shows the power spectrum from the simulations for model M001 at z = 0. The data have three notable features that are important for the modeling procedure. (1) The non-standard representation for the power spectrum

$$\begin{equation} {\cal P}(k)=\Delta ^2(k) /k^{1.5} \end{equation} \tag{ 5 }$$

is chosen to accentuate the BAOs. We will account for this feature when we choose our function class for the smooth power spectra (discussed later). (2) The three series of simulations G, H, and L with 1, 4, and 16 realizations, respectively. Because the H- and L-series do not have enough force resolution to resolve the nonlinear regime at high k, we need to restrict the use of these runs to small and intermediate k ranges. (3) The simulation variance at any given k is known.

Zoom In Zoom Out Reset image size

We treat each simulated ("noisy") realization of the (large volume or averaged) spectrum as a draw from a multivariate Gaussian distribution whose mean is given by an unknown smooth spectrum. Thus, for a given cosmology c, a given series s ∈ {G, H, L}, and a given replicate i = 1, ..., N_s (where N_s is the number of simulations for the series), we have a multivariate Gaussian density for the simulated spectrum P^c_s,i,

$$\begin{eqnarray} f\big(P^{c}_{s,i}\big) & \propto & \left| A_s \Omega A_s^{T} \right|^{1/2}\nonumber\\ && \times \exp \bigg\lbrace -\frac{1}{2} \left(P^{c}_{s,i} - A_s \mathcal {P}^{c} \right)^{T} \nonumber\\ && \times A_s \Omega A_s^{T} \left(P^{c}_{s,i} - A_s \mathcal {P}^{c} \right) \bigg\rbrace. \end{eqnarray} \tag{ 6 }$$

Here, $\mathcal {P}^c$ is the smooth power spectrum; A_s is a projection matrix of zeros and ones that is used to remove the high-k values for which the H- and L-series are not used (thus, A_G is an identity matrix); and Ω is a diagonal matrix of the known precisions (inverse variances).

The model is completed by specifying a class of functions for the smooth power spectrum $\mathcal {P}^c$. We choose a flexible class of smooth functions called a process convolution (Higdon 2002). These functions are best described constructively. A process convolution builds a smooth function as a moving average of a simple stochastic process like independent and identically distributed (i.i.d.) Gaussian variates or Brownian motion. The moving average uses a smoothing kernel whose width is allowed to vary over the domain to account for nonstationarity (smoother in some regions, more wiggly in others). Figure 3 shows a simple example of a process convolution built on white noise smoothed with a Gaussian kernel.

Zoom In Zoom Out Reset image size

We build the process convolution for $\mathcal {P}^c$ on Brownian motion u^c, with marginal variance τ²_u, realized on a sparse grid (relative to the power spectrum) of evenly spaced points, x (the number of points is not important so long as it is large enough; we use 100),

$$\begin{equation} f(u^c) \propto \left| \frac{1}{\tau ^2_u} W \right|^{1/2} \exp \left\lbrace -\frac{u^{c \prime } W u^c}{2 \tau ^2_u} \right\rbrace, \end{equation} \tag{ 7 }$$

where W is the Brownian precision matrix with diagonal equal to [12...21] and −1 on the first off-diagonals (this matrix cannot actually be inverted, but never has to be in the estimation). The Brownian motion is transformed into $\mathcal {P}^c$ by the smoothing matrix K^σ,

$$\begin{equation} \mathcal {P}^c = K^{\sigma } u^c. \end{equation} \tag{ 8 }$$

The smoothing matrix K^σ is built using Gaussian kernels whose width varies smoothly across the domain. Thus, we have

$$\begin{equation} K^{\sigma }_{i,j} = \frac{1}{\sqrt{2 \pi \sigma ^{2}_{i}}} \exp \left\lbrace -\frac{ \left(\log _{10}(k_i) - x_j \right)^2}{2 \sigma _i^2} \right\rbrace, \end{equation} \tag{ 9 }$$

where k_i is the ith value of k for which the power spectrum is computed.

In the description of K^σ, σ is indexed to indicate that it changes over the domain. Intuitively, we want σ to be small in the middle of the domain in order to capture the oscillations, but large elsewhere to smooth away the noise. In order to estimate this varying bandwidth parameter, we build a second process convolution model. This model is built on i.i.d. Gaussian variates v, with mean zero and variance τ²_v, observed on an even sparser grid of evenly spaced points, t (length M_v; we use 10, but any large enough number will suffice),

$$\begin{equation} f(v) \propto \left(\frac{1}{\sqrt{ \tau _{v}^{2} } } \right)^{M_v} \exp \left\lbrace -\frac{v^{\prime } v}{2 \tau ^2_v} \right\rbrace. \end{equation} \tag{ 10 }$$

The process v is transformed into σ by the smoothing matrix K^δ,

$$\begin{equation} \sigma = K^{\delta } v. \end{equation} \tag{ 11 }$$

This matrix is also built using Gaussian smoothing kernels, but with a constant bandwidth, δ. Thus, we have

$$\begin{equation} K^{\delta }_{i,j} = \frac{1}{\sqrt{2 \pi \delta ^{2}}} \exp \left\lbrace -\frac{\left(x_i - t_j \right)^2}{2 \delta ^2} \right\rbrace. \end{equation} \tag{ 12 }$$

Combining all of this, we get a distribution for the simulated power spectra for a given cosmology,

$$\begin{eqnarray} f\big(P^{c}_{s,i}\big) & \propto & \left| A_s \Omega A_s^{\prime } \right|^{1/2}\nonumber\\ && \times \exp \bigg\lbrace -\frac{1}{2} \left(P^{c}_{s,i} - A_s K^{\sigma } u^c \right)^{\prime } \nonumber\\ && \times A_s \Omega A_s^{\prime } \left(P^{c}_{s,i} - A_s K^{\sigma } u^c \right) \bigg\rbrace, \end{eqnarray} \tag{ 13 }$$

where s ∈ {G, H, L} and i = 1, ⋅ ⋅ ⋅, N_s. Note that K^σ is a function of the parameters v and δ. We choose noninformative priors for τ²_u, τ²_v, and δ so as to impart little or no information about their values,

$$\begin{eqnarray} \pi \big(\tau ^2_u \big) & \propto & \left(\frac{1}{\tau ^2_u} \right)^2 \exp \left(-\frac{0.001}{\tau ^2_u} \right), \nonumber \\ \pi \big(\tau ^2_v \big) & \propto & \left(\frac{1}{\tau ^2_v} \right)^2 \exp \left(-\frac{0.001}{\tau ^2_v} \right),\nonumber\\ \pi (\delta) & \propto & I\ \ \lbrace 0 \le \delta \le 10 \rbrace, \end{eqnarray} \tag{ 14 }$$

which are inverse gamma, inverse gamma, and uniform, respectively.

Equations (7) (for all c), (10), and (14) are multiplied together with Equation (13) (for all c) to produce a posterior distribution for the unknown parameters and stochastic processes. Obtaining an estimate for the smooth power spectra requires an estimate for each of the parameters, the process v, and each of the processes u^c for all c. Markov Chain Monte Carlo (MCMC) via the Metroplis–Hastings algorithm (Chib & Greenberg 1995) produces a sample from the posterior distribution by drawing each parameter individually. The stochastic processes u^c can be integrated out of the distribution, so the MCMC produces samples of the parameters as well as v. We use the posterior mean of the parameters to obtain a conditional mean for each of the u^c which is then transformed into each of the $\mathcal {P}^c$.

The perturbation results are included by setting the low-k values for each of the simulated spectra in a given cosmology to the perturbation results and setting the precisions for these values to be very large. This is done prior to the estimation described here. The replication of these values in every simulation, combined with the large precisions, nearly forces the estimated result through these points. Further, the transition from the N-body results to the perturbation results will be relatively smooth as long as the two line up fairly well.

3.2.2. Tests with the Linear Power Spectrum

The next step is to test the matching procedure between theoretical and simulation results to ensure high accuracy at the perturbation theory matching point. Since we do not know the exact answer for the nonlinear power spectrum, we first carry out a test using linear theory. In this test, we use the power spectra from the initial conditions and show how well we can smooth out the run-to-run scatter. Knowing the exact answer here will allow us to assess how well the matching procedure actually works.

For these spectra, the entire estimation procedure is applied to the initial condition spectra. The only difference is that fewer grid points were used for the latent processes u and v to account for the reduced range of the spectra (we use 70 and 7, respectively). The final results are shown in Figure 4 for the same set of three random models as in Figure 1. We verified that the results hold for all of the remaining models. The upper panels of each sub-panel show the theoretical linear power spectrum in black and the prediction from the simulations in red. The vertical line marks the matching point to linear theory at k ≃ 0.03 Mpc⁻¹. The lower panels show the ratio of the predicted power spectra to the theoretical power spectra in red. Below the matching point, the agreement is—by construction—perfect. Beyond the matching point, the smoothed prediction in red is accurate at the 1% level. This test shows that we can obtain a smooth, high-accuracy prediction for the power spectrum by combining a suite of realizations with perturbation theory.

Zoom In Zoom Out Reset image size

3.3. The Nonlinear Power Spectrum

The process convolution procedure is applied to the power spectrum realizations for each of the 37 cosmologies at six values of the scale factor a = 1/(1 + z) ∈ {0.5, 0.6, 0.7, 0.8, 0.9, 1.0}, where z is the redshift. The results are shown in Figure 5, again for a subset of three models. Although we have no known truth to use as a comparison in this case, the resulting estimates continue to fit the simulation realizations very well.

Zoom In Zoom Out Reset image size

In place of comparing with known truth, we can examine some tests of our modeling assumptions. We assume that the simulation spectra on the modified scale are independently and normally distributed about a smooth mean with a variance that changes with k. Given this, we can compute standardized residuals in which the estimated smooth mean is subtracted from the simulations and the result is scaled by the known standard deviation. These standardized residuals should look like i.i.d. standard normal variables. Figure 6 shows quantile–quantile plots for the G simulations of three cosmologies at each of the six scale factors. Theoretical quantiles from a standard normal distribution are given on the x-axis and the sample quantiles of the standardized residuals are shown on the y-axis. The nearly straight line at 45° indicates an extremely good distributional fit. Figure 7 shows the standardized residuals plotted against k for the same simulations. This plot verifies the independence assumption. Further, there is an almost complete lack of evidence of any structure in these plots, suggesting that we would not greatly improve the fit by relaxing the smoothness assumption.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is also interesting to examine the plot of the kernel width function as estimated by the MCMC process. Figure 8 shows the median draw for σ as a function of k. As expected, the kernel width is small in the vicinity of the baryon wiggles. This means that the values of the latent process u in this region have large weights relative to values further away. This prevents real local structure, like the BAO, from being smoothed out. On both ends, the kernel width is large, which results in a smooth function in these regions.

Zoom In Zoom Out Reset image size

Having extracted the smooth power spectra from our simulation suite, we can now build an emulator to predict the nonlinear matter power spectrum within the priors specified in Equations (2). We will use only models 1–37 for the emulator construction; model M000 will serve as an independent check of the emulator accuracy, along with hold-out tests described below.

In order to construct the emulator, we model the 37 power spectra using an $n_\mathcal {P}$ dimensional basis representation:

$$\begin{equation} \mathcal {P}(k;z;\theta)= \sum _{i=1}^{n_\mathcal {P}}\phi _i(k;z)w_i(\theta), \theta \in [0, 1]^{n_\mathcal {\theta }}, \end{equation} \tag{ 15 }$$

where the ϕ_i(k; z) are the basis functions, the w_i(θ) are the corresponding weights, and the θ represent the cosmological parameters. The dimensionality $n_\mathcal {P}$ refers to the number of orthogonal basis vectors $\lbrace \phi _i(k,z),\dots,\phi _{n_\mathcal {P}}(k,z)\rbrace$. The parameter n_θ is the dimensionality of our parameter space—with five cosmological parameters we have n_θ = 5. The power spectrum $\mathcal {P}(k;z;\theta)$ depends on the wavenumber k, the redshift z, and the five cosmological input parameters θ (note that we rescale the range of each parameter to [0, 1]). Examination of the results indicates that $n_\mathcal {P}=5$ is a good choice for the number of basis vectors ϕ_i(k; z) (that $n_\mathcal {P}=n_\theta$ here is a coincidence). The task is now to (1) construct a suitable set of orthogonal basis vectors ϕ_i(k; z) and (2) model the weights w_i(θ)). For the first task we use principal components, and for the second, Gaussian Process (GP) models. Our choice of GP modeling is based on their success in representing functions that change smoothly with parameter variation, e.g., the variation of the power spectrum as a function of cosmological parameters. Both steps are explained in great detail in Paper II; we refer the interested reader to that publication. Paper II also contains various error control tests of the GP-based interpolation method which we do not repeat here.

Since the details of how to build an emulator are already provided in Paper II we can immediately turn to our final product: the emulator itself. To facilitate use of the emulator, we are releasing the fully trained emulator with this paper. It provides nonlinear power spectra at a set of redshifts between z = 0 and z = 1 out to k = 1 h Mpc⁻¹, for any cosmological model specified within the priors given in Equations (2).

4.1. Emulator Test

In order to verify the accuracy of the emulator, we perform two important checks (other tests on the methodology were carried out in Paper II). For the first, we compare the simulation results for model M000 with the emulator prediction. Figure 9 shows the results of this true out-of-sample prediction. The plot shows the ratio of the emulated power spectrum to the actual simulation for the M000 cosmology at six values of the scale factor a. This cosmology is completely interior in the design and, since M000 was not used to estimate the smoothing or the emulator, this test provides a good indication of the actual performance of the emulator within its parameter bounds. Overall the agreement between the simulations and the emulator is excellent and for most of the k-range well below the 1% error bound.

Zoom In Zoom Out Reset image size

The second check consists of a sequence of holdout tests. In a holdout test, the emulator is built from 36 cosmological models and the emulator prediction can be then compared to the result from the extra model. The drawback of this test is that if we have only a very small number of models, each of them is important for capturing some part of the parameter space and taking it out for building the emulator degrades the emulator precision. In order to keep this problem to a minimum, we only perform holdout tests for models which are interior simulations, meaning that none of the five parameters are near the extreme limits of the chosen prior range. There are six such models: M004, M008, M013, M016, M020, and M026. The ratio of the emulated prediction to the actual simulation for these models is shown in Figure 10. For each cosmology, there are six ratios for each value of the scale factor. The lines are quite well behaved with errors largely within the 1% bounds.

Zoom In Zoom Out Reset image size

The advent of precision cosmology and the prospect of very large surveys such as LSST and JDEM poses an enormous challenge to the theory community in the field of large-scale structure predictions. Future progress not only requires very accurate predictions but also the ability to produce predictions for different cosmological models very fast. The aim of the Coyote Universe project was to take a first step in attacking this problem, focusing on the matter power spectrum at intermediate scales out to k ≃ 1 h Mpc⁻¹. While predicting the matter power spectrum on these scales at high accuracy may appear to be a moderately difficult task, it turned out to be a technical and computational challenge in several respects. It is generally hard to imagine problems and pitfalls in advance, which is why fully working through an example is so helpful. Since the community has to follow a similar path in the future to create predictive capabilities for different cosmological probes, we summarize here some of the lessons learned during this project.

The major differences between a project like this (including all three papers of the series) and previous numerical studies of large-scale structure probes are:

• 1.  
  Computational and Storage Capacity. The Coyote Universe simulation suite encompasses roughly 60 TB of data. The computational cost is of the order of a million CPU-hr and including waiting times in submission queues, downtimes of the machine, and so on, carrying out these simulations took roughly six months. The simulation size of one billion particles in a Gpc³ volume was barely enough to resolve the scales of interest and for future work will certainly not be sufficient. Such simulations will need larger volumes and better mass resolution. For example, to resolve a 10¹² M_☉ halo with 100 particles in a (3 Gpc)³ volume, we would need 300 billion particles. While supercomputers will get faster and larger in the future, generating many simulations at the edge of machine capabilities will always be a challenge. In addition, archiving the outputs of such simulations will become very expensive in terms of storage. From the Coyote Universe runs we stored 11 time snapshots plus the initial conditions (particle positions and velocities and halo information) leading to 250 GB of data per run. For the 300 billion particle run this would increase to 75 TB. Only very few places worldwide would be able to manage the resulting large databases.
• 2.  
  Simulation Infrastructure. Running a very large number of simulations makes it necessary to integrate the major parts of the analysis steps into the simulation code and to automate as much of the mechanics of running the code (submission, restarts) as possible. For the Coyote Universe project we developed several scripts to generate the input files of the codes, to structure the directories where different runs are performed in, and to submit the simulations to the computing queue system. For future efforts of this kind the adoption and development of dedicated workflow capabilities for these tasks must be considered. The number of tasks to be carried out will become too large to keep track of without such tools. In addition, since large projects will require extensive collaborations, software tools will make it easier to work in a team environment since each collaborator will have information about previous tasks and results. An example of such a tool for cosmological simulation analysis and visualization is given in Anderson et al. (2008).We carried out the data analysis after the runs were finished. For very large simulations this is not very practical, and on-the-fly analysis tools are required to minimize read and write times and failures. This in turn requires that the code infrastructure be tailored to the problem under consideration.
• 3.  
  Serving the Data. Clearly, large simulation efforts cannot be carried out by a few individuals, and require possibly community-wide coordination. The simulation data will be valuable for many different projects. It is therefore necessary to make the data from such simulation efforts publicly available and serve them in a way that new science can be extracted from different groups of simulations. Transferring large amounts of data is difficult because of limitations in communication bandwidth and also because of the large storage requirements. It would therefore be desirable to have computational resources dedicated to the database. In such a situation, researchers would be able to run their analysis codes on machines with direct access to the database and perform queries on the data easily. We are planning to make the Coyote Universe database available in the future and use it as a manageable testbed for such services.
• 4.  
  Communication with Other Communities. The complexity of the analysis task makes it necessary to efficiently collaborate and communicate with other communities, for example, statisticians, computer scientists, and applied mathematicians. Many tools that will be essential for precision cosmology in the future have already been invented—the task is to find them and use them in the best way possible.

This paper is the last of the Coyote Universe series of publications. Paper I was concerned with demonstrating that percent level accuracy in the (gravity only) nonlinear power spectrum could be attained out to k ≃ 1 h Mpc⁻¹. Paper II showed that with only a relatively small number of simulations, interpolation across a high-dimensional space was possible at close to the same level of accuracy as that attained in the individual runs. Paper III takes this work to the final conclusion: based on almost 1000 simulations spanning 38 wCDM cosmologies, we present a fast and very accurate prediction scheme—an emulator—for the nonlinear matter power spectrum. The emulator is accurate at the percent level, improving over commonly used fitting functions by almost an order of magnitude.

The emulator construction—as explained in Paper II—is based on GP modeling. In order to carry this out, a major challenge is to produce a smooth power spectrum from a finite set of simulations for each cosmology. In order to minimize run-to-run scatter on very large scales, we performed several medium-resolution simulations and matched these at sufficiently low k to perturbation theory. On quasi-linear scales we used medium-resolution simulations and matched those to high-resolution runs at small spatial scales. Matching the different resolution runs and perturbation theory accurately was carried out using process convolution. This technique allowed us to construct a smooth power spectrum for each cosmological model. The results were then used to construct the emulator via GP modeling as described in detail in Paper II.

We are releasing the power spectrum emulator as a C code, which allows the user to specify a cosmology within our priors and returns the power spectrum at six different redshifts between z = 0 and z = 1 out to k ≃ 1 h Mpc⁻¹. These power spectra can now be used for further analysis of cosmological data. We are planning to extend the emulator in the near future to a larger k-range and will provide a smooth interpolation between results from different redshifts.

A major challenge will be to ensure that a certain level of accuracy is reached with the simulations when going to small scales. Besides being computationally very expensive (high force and mass resolution being required) the physics at smaller scales is far more complicated. The inclusion of gasdynamics and feedback effects (along with other physics) is far from being straightforward. The impossibility of carrying out a direct simulation effort is more or less certain; a good number of phenomenological/subgrid modeling parameters will be required. As simulation complexity and the number of modeling and cosmological parameters increases, it becomes even more important to develop efficient and controlled sampling schemes as described in Paper II, so that data can be used to determine both cosmological and modeling parameters (self-calibration).

With the series of three Coyote Universe papers we have demonstrated that it is possible to extract cosmological statistics such as the power spectrum at high accuracy and that one can build an accurate prediction scheme based on a limited set of simulations. This line of work will be important for interpreting results of future cosmological surveys. It will also have to be extended in several ways: (1) the cosmological model space has to be opened up; (2) we have to ensure high accuracy at scales smaller than those considered here; (3) we have to include more physics in order to capture those small scales correctly; and (4) we have to include different cosmological probes, e.g., the cluster mass function and the shear power spectrum, to be able to build a complete framework for analyzing future survey data. We have shown here that such a program can in principle be established though it will demand a large concerted effort between different communities.

A special acknowledgment is due to supercomputing time awarded to us under the LANL Institutional Computing Initiative. Part of this research was supported by the DOE under contract W-7405-ENG-36 and by a DOE HEP Dark Energy R&D award. S.H., K.H., D.H., E.L., and C.W. acknowledge support from the LDRD program at Los Alamos National Laboratory. K.H. was supported in part by NASA. M.W. was supported in part by NASA and the DOE. We thank Dragan Huterer, Nikhil Padmanabhan, Adrian Pope, and Michael Schneider for useful discussions. We thank Volker Springel for making the N-body code GADGET-2 publicly available.