Dark Quest. I. Fast and Accurate Emulation of Halo Clustering Statistics and Its Application to Galaxy Clustering

In Figure 5, we first examine the HMF for the fiducial Planck cosmology at three redshifts, z = 1.02, 0.549, and 0, using four N-body simulations with different numerical resolutions. Note that the simulation with the lowest resolution among the four, which has 256³ particles, has the same resolution as our main LR suite, whereas the second from the worst, with 512³ particles, corresponds to the resolution of the HR suite. For reference, the HMF in Figure 5 is normalized by the fitting formula of Tinker et al. (2008), with a mass definition of 200 times the cosmic mean density. To have a fair comparison, we integrated the Tinker et al. (2008) HMF (hereafter Tinker HMF) over the halo masses in each mass bin, which are used when we measure the HMF from simulations.

Zoom In Zoom Out Reset image size

We can see that the measured HMF better matches the Tinker HMF down to lower masses as the simulation resolution increases. The four simulations agree with each other at the high-mass end, although the curves are noisy due to the Poisson noise. Note that these simulations are done in a comoving box with a side length of 250 h⁻¹ Mpc, which is smaller than our main simulations of 1 or 2 h⁻¹ Gpc used for the emulator development. This suggests that our main simulations have much lower Poisson noise at such high-mass bins. The vertical arrows, from right to left for higher resolution, denote the halo mass corresponding to 100 N-body particles. The figure indicates that the simulation HMF at this mass scale is underestimated by about 10%, fairly independently of redshift. Thus, one needs at least several hundred particles to estimate the HMF to a percent accuracy.

We here propose a way to empirically correct for a systematic error in the estimated HMF due to numerical resolution. Our method is motivated by the method in Warren et al. (2006), which was developed for halos that are identified by the friends-of-friends (FoF) method. They proposed that the FoF mass of each halo is calibrated as

$$\begin{eqnarray}&&\tilde{M}=\left(1-{N}_{{\rm{p}}}^{-0.6}\right)M,\end{eqnarray} \tag{ 26 }$$

where N_p = M/m_p is the number of member particles, m_p is the N-body particle mass, and $\tilde{M}$ is the corrected mass. Since the FoF algorithm tends to link physically unbound particles near the halo boundary when the mass resolution is poor, an FoF halo mass tends to be overestimated compared to the true mass. Hence, the FoF-based HMF tends to be overestimated for low halo masses that are affected by numerical resolution. The correction factor is applied to each FoF halo in such a way that the FoF mass is reduced to correct for the overestimation in HMF. This procedure was further confirmed in Crocce et al. (2010), where the method was applied to FoF halos in MICE simulations. On the other hand, our result in Figure 5 displays a rather opposite trend: our HMF is underestimated in low-resolution simulations, implying that the mass of each low-mass halo tends to be underestimated. Since we use the SO mass, the SO-based HMF is affected by the matter density field around the halo region, which tends to be underestimated in low-resolution simulations. This is different from the FoF finder, and thus the opposite trend is understandable.

We thus use the following equation to correct for the SO mass:

$$\begin{eqnarray}&&\tilde{M}=\left(1+{N}_{{\rm{p}}}^{-0.55}\right)M,\end{eqnarray} \tag{ 27 }$$

where N_p is the number of particles within R_200m in the SO mass definition. We employ a slightly different power of N_p from Equation (26) to get a better calibration. After this correction, the HMFs from different resolution simulations better agree with each other down to the resolution limit denoted by the vertical arrows, as shown in Figure 6. Below the mass limit, the SO halo mass becomes overcorrected, yielding an overestimation in HMF. These trends after the correction appear to be very similar at different redshifts. While a further refinement of the empirical function given by Equation (27) would be possible in principle, we do not use halos with less than 200 particles when we calibrate HMF. Note that we employed a slightly conservative threshold of 200 particles, as compared to the case of 100 particles discussed in Figure 6. Furthermore, we use only the HR simulations (resolution equivalent to the one with 512³ particles here) for development of the HMF emulator. We can determine HMF accurately down to ∼10¹² h⁻¹ M_⊙, and this number varies depending on the cosmological model, as we will discuss below.

Zoom In Zoom Out Reset image size

Next, we check the clustering correlation functions of halos. In doing so, we need to consider subsamples of halos divided by halo discriminators, such as halo mass, and then consider the clustering correlations as a function of different subsamples. In the left plot of each panel in Figures 7 and 8, we show the halo–matter cross- and halo–halo autocorrelation functions for a mass threshold sample of halos with $M\geqslant {10}^{12}\,{h}^{-1}$ or $5\times {10}^{12}\,{h}^{-1}{M}_{\odot }$ measured from the four simulations as in Figure 5. Here the threshold mass ${10}^{12}\,{h}^{-1}{M}_{\odot }$ corresponds to halos with more than 100 member particles for simulations with more than 512³ particles, whereas it corresponds to halos with only ∼10 particles for the 256³ simulation. Note that we did not apply the correction from Equation (27) for halo masses in these figures. While all of the correlation functions agree with each other at large separations, the smaller scales clearly show the effect of numerical resolution; the measurements from a lower-resolution simulation start to deviate from those from a higher-resolution simulation on scales smaller than ∼1 h⁻¹ Mpc. The comparison of Figures 7 and 8 reveals that the deviation is larger for a halo sample of smaller mass threshold. The inaccuracy is ascribed to several facts. In a lower-resolution simulation, halo masses around the mass threshold are not determined accurately on an individual halo basis due to the lack of numerical resolution, as discussed in Figure 5. Thus, the halo sample of a given mass threshold becomes different from that of a higher-resolution simulation. Moreover, the mass distribution around each halo in a lower-resolution simulation is simulated less accurately.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this paper, we employ a slightly different sample of halos to develop the emulator. Rather than using the mass as the primary proxy of the different clustering strengths of halos, we consider mass threshold samples and label each sample in terms of the number density of halos above the threshold. We expect two advantages from this conversion. First, the cosmology dependence of the noise level of various statistics is weaker compared to the samples labeled by the mass. Indeed, we know that the mass of the heaviest halos available in each simulation can be quite different among different cosmological models and at different redshifts. Second, as is quite obvious from Figure 5, the masses inferred from simulations are quite sensitive to the resolution, especially at the low-mass end. To see this more qualitatively, we show in the right plot of each panel in Figures 7 and 8 the clustering signals for the mass threshold samples with a fixed number density in each simulation. Here the number density is the same as that of the mass threshold sample for the highest-resolution simulation (2048³) in the left plot, and the mass thresholds for other simulations are determined to match the target number density. Now an agreement between different resolution simulations is better than in the left panel, reflecting the fact that the number of halos in the sample is less affected by numerical resolution compared to the mass threshold sample. Nevertheless, the lowest-resolution simulation still exhibits a relatively large deviation for the halo–matter cross-correlation at small scales, especially for the sample corresponding to ${10}^{12}\,{h}^{-1}{M}_{\odot }$, because the matter distribution in high-density regions is less accurately simulated in such a low-resolution simulation. To avoid this inaccuracy, we use only the HR simulations to estimate the halo–matter cross-correlation function for different cosmological models.

On the other hand, Figure 8 shows a slightly better agreement among the four simulations with different resolutions, implying that the halo–halo autocorrelation is relatively robust against the numerical resolution. Note that larger scatters in the ratio on small scales ($\lesssim 1\,{h}^{-1}\mathrm{Mpc}$) are due to the halo exclusion effect, which states that the correlation signal is sharply suppressed on small scales due to the fact that no halo pair can exist below R₂₀₀ radius of the larger one by construction in our halo sample. Thus, a slight misestimation of the halo radii due to numerical resolution can lead to a large error in the correlation signal at a fixed scale around the typical R₂₀₀ of the sample. Since our final product is the galaxy correlation function, and the one-halo term gives a dominant contribution around these scales, the scatter seen here does not largely affect the predictions of the galaxy autocorrelation function, as we will show later. Based on these results, we use the LR simulations to estimate the autocorrelation functions for different cosmological models.

In summary, we use the correlation functions of halos measured for halo samples with different number densities. When this is combined with the HMF module that gives the halo number density as a function of mass, one can compute the halo correlation function for a given mass threshold instead of the number density. Furthermore, we can compute the correlation function of halos in an infinitesimally narrow mass bin by taking the numerical derivative of the correlation functions for a mass threshold halo sample with respect to the threshold mass.

The clustering correlation functions of halos measured from each of the simulations become considerably noisy on very large scales around the baryon acoustic oscillation (BAO) scale due to the significant sample variance due to the finite simulation volume, even for our LR simulations of 2 h⁻¹ Gpc size. To overcome this obstacle, we employ a semi-analytical approach based on the propagator (e.g., Crocce & Scoccimarro 2006a), which captures most of the expected linear and nonlinear effects around the BAO scale. We then stitch the prediction with the direct simulation results to obtain model predictions over a wide range of scales, as described below.

Figures 9–11 show the matter auto-, halo–matter cross-, and halo autocorrelation functions on large scales, respectively. The solid curves in the upper left panel of each figure depict the correlation functions measured from each of the 14 LR realizations for the fiducial Planck cosmology. Clearly, the realization-to-realization scatter is large. For comparison, the upper right panels show the linear theory predictions we computed using the same Gaussian random realizations as in the initial conditions of each simulation in order to properly take into account the sample variance effect (see Section 4.3.4 for our method to determine the linear bias parameter from the halo correlation functions). The scatter among the realizations seen in the linear predictions is comparable to that of the corresponding nonlinear counterparts, except for the halo–halo autocorrelation function with a low number density of ${10}^{-5}\,{({h}^{-1}\mathrm{Mpc})}^{-3}$ (i.e., the right panels of Figure 11). This suggests that the primary source of the scatter is indeed the sample variance in the initial conditions, and the shot noise adds only a moderate scatter for low-density samples of halos.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since our varied cosmology simulation suite is, in principle, performed only once at each model, the large scatters in the measured correlation function make it difficult to construct an accurate emulator. Unlike the matter power spectrum, we cannot switch to a parameter-free perturbative calculation on these large scales because we have to know the halo bias that is not accurately described by a simple analytical prescription that often ignores the dependence on scale and cosmology. We need an appropriate method where the sample variance is sufficiently reduced and, at the same time, the large-scale bias of the halos under consideration is properly taken into account.

Another important effect on the BAO scale, in addition to the large-scale bias, is the damping of the BAO feature. This is clearly visible from comparison of the upper left and right panels in Figures 9–11. It is known that this effect is to a large extent due to the large-scale bulk motion of the cosmic fluid, which can be accurately modeled by the propagator (Crocce & Scoccimarro 2006a). In their paper, the propagator for the matter field is defined by

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{\partial {{\rm{\Phi }}}_{a,{\boldsymbol{k}}}}{\partial {\delta }_{\mathrm{lin},{\boldsymbol{k}}^{\prime} }}\right\rangle \equiv {\delta }_{{\rm{D}}}^{3}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ){G}_{a}(k),\end{eqnarray} \tag{ 28 }$$

where Φ_a can be either the density or the velocity divergence field of matter. Note here and in what follows that the linear density field δ_lin and its power spectrum P_lin are always scaled by the linear growth factor to the same redshift as other quantities, such as Φ_a or G_a. The function G_a(k) is called the (two-point) propagator, which shows a damping form very close to a Gaussian shape toward high k. This function can be interpreted to describe how much memory of the initial density field (δ_lin) persists in the final (nonlinear) fields (Φ_a). One can analytically show that this function is exactly a Gaussian with its variance equal to the inverse square of the rms displacement field in the case of the Zel'dovich dynamics for a Gaussian initial condition.

In most of the resummed perturbation theories, the leading-order contribution to the mixed power spectrum of two fields δ_a and δ_b is expressed as ${G}_{a}(k){G}_{b}(k){P}_{\mathrm{lin}}(k)$, where the subscripts a and b can be the density or velocity divergence of matter or any tracers (e.g., Crocce & Scoccimarro 2006b, 2008; Bernardeau et al. 2008) and P_lin(k) is the linear matter power spectrum. The IFT of this combination gives a reasonable prescription on the two-point correlation function around the BAO scale,

$$\begin{eqnarray}&&{\xi }_{a,b,\mathrm{tree}}(r)=\mathrm{iFT}\left[{G}_{a}(k){G}_{b}(k){P}_{\mathrm{lin}}(k)\right],\end{eqnarray} \tag{ 29 }$$

where we use the subscript "tree" to indicate that this quantity is the tree-level result (i.e., the leading-order diagrams) of the resummed perturbation theories. Indeed, Equation (29) with a simple Gaussian approximation of the propagator can already explain the damping of the BAO peak in the matter correlation function very accurately (e.g., Matsubara 2008). By going into higher orders, a subpercent-level shift in the BAO peak location to a smaller separation scale can be realized (Crocce & Scoccimarro 2008). This will be important in interpreting the BAO-related distance measurements from actual observations.

We now consider the propagator for halos. One can define the propagator by simply replacing Φ_a with the density field of halos in a given sample. In what follows, we denote this by G_a(k), where the subscript a is either the matter or the halo density field. In case of halos, the low-k limit of the function corresponds to the linear bias factor. A damping behavior at high k should be similar to that of the matter field, and this damping is responsible for the smearing of the BAO peak measured through the clustering of halos. We show the functions for matter and halos with different number densities and at different redshifts in Figure 12, and they are measured from the 14 LR realizations of the fiducial Planck cosmology. We estimate the function by taking

$$\begin{eqnarray}&&{G}_{a}(k)=\displaystyle \frac{{P}_{a,\mathrm{lin}}(k)}{{P}_{\mathrm{lin}}(k)},\end{eqnarray} \tag{ 30 }$$

where ${P}_{a,\mathrm{lin}}(k)$ is the cross-power spectrum of tracers "a" and the linear density field. In this estimator, we use the linear power spectrum P_lin(k) measured from the linear density field used for the initial condition instead of the theoretical smooth function, and the sample variance is largely suppressed by taking this ratio. Indeed, the scatter among the 14 realizations seen in the figure is rather small, especially for low number density halo samples compared to the scatter in the corresponding autocorrelation function in Figure 11. The overall trend of this function looks very simple, as already discussed; it appears to be a Gaussian-like damping function with a linear bias factor at the low-k limit that depends on the halo number density. In addition, we can see that the damping starts at smaller k at lower redshifts, reflecting the fact that the information in the initial density field remains more on larger scales and at higher redshifts.

Zoom In Zoom Out Reset image size

To summarize, our strategy to describe the large-scale limit of matter or halo correlation functions is to emulate the function G_a(k) for both matter and halo fields and substitute it into Equation (29). Likewise, we take the same combination for the random fields δ_lin used in the initial conditions, which we schematically denote as $\mathrm{iFT}[{(G{\delta }_{\mathrm{lin}})}^{2}]$. We show this model in the lower left panel of Figures 9–11 for the random realizations corresponding to the 14 simulations shown in the upper panels. The curves obtained in this way appear to be very similar to the direct simulation results in the upper left panels.

Finally, the averages of these curves are shown by the downward-pointing triangles with error bars in the lower right panel. They are almost indistinguishable from the circles for the nonlinear correlation function directly measured from the nonlinear fields. Indeed, their differences, shown by the crosses, are consistent with zero. A closer look at the scale dependence of this residual indicates a small pattern that would cause a small shift on the BAO peak toward a smaller scale. It tends to be positive around the inflection point of the correlation function (at around 90 h⁻¹ Mpc) and negative at scales smaller than 80 h⁻¹ Mpc. Since in most cases, these features are within the error bars, which correspond to the scatter among realizations, we simply ignore this small residual in the following discussion.

We also show the continuous limit of the model, Equation (29), by the dashed line. This is the expectation value of the downward-pointing triangles in the limit of an infinite number of realizations. Our final model for the large-scale correlation function is this line. With this procedure, we can reduce the sample variance significantly, since the prediction is based on the noiseless linear power spectrum P_lin. Our approach works well even in the case of the halo autocorrelation function for a halo sample with a small number density (see the upper left panel in the right part of Figure 11); the unaccounted shot-noise effect only adds a random scatter, and no systematic trend can be seen in the residual. We will explain how we switch from the direct measurement of the correlation functions to the prescription based on the propagator explained here in Section 5.1.1.

In this section, we describe how to build a module of the HMF. As shown in Figures 5 and 6, the fitting formula by Tinker et al. (2008) works very well, at least for the fiducial cosmology at z = 0. The fitting function we use in the following is a modified version of the earlier model in Press & Schechter (1974; see also Sheth & Tormen 2002), given by

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{dM}}=f({\sigma }_{M})\displaystyle \frac{{\bar{\rho }}_{{\rm{m}}}}{M}\displaystyle \frac{d\mathrm{ln}{\sigma }_{M}^{-1}}{{dM}},\end{eqnarray} \tag{ 31 }$$

with

$$\begin{eqnarray}&&f({\sigma }_{M})=A\left[{\left(\displaystyle \frac{{\sigma }_{M}}{b}\right)}^{a}+1\right]\exp \left(-\displaystyle \frac{c}{{\sigma }_{M}^{2}}\right).\end{eqnarray} \tag{ 32 }$$

Here the mass variance ${\sigma }_{M}^{2}$ is given by

$$\begin{eqnarray}&&{\sigma }_{M}^{2}=\int \displaystyle \frac{{k}^{2}{dk}}{2{\pi }^{2}}{P}_{\mathrm{lin}}(k;z){\left|{\tilde{W}}_{R}(k)\right|}^{2},\end{eqnarray} \tag{ 33 }$$

where P_lin(k; z) is the linear matter power spectrum at redshift z, and ${\tilde{W}}_{R}(k)$ is the Fourier transform of a top-hat filter of radius R that is specified by an input halo mass M via $R\,={(3M/4\pi {\bar{\rho }}_{{\rm{m}},0})}^{1/3}$. Tinker et al. (2008) showed that the HMF measured in the simulations is well fitted by the above functional form with time-dependent coefficients:

$$\begin{eqnarray}&&A(z)=0.186{\left(1+z\right)}^{-0.14},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&a(z)=1.47{\left(1+z\right)}^{-0.06},\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&b(z)=2.57{\left(1+z\right)}^{-\alpha },\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&c(z)=1.19,\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&\alpha =-{\left(\displaystyle \frac{0.75}{{\mathrm{log}}_{10}({\rm{\Delta }}/75)}\right)}^{1.2}.\end{eqnarray} \tag{ 38 }$$

The overdensity Δ is 200 in our halo mass definition.

The variation in the HMF over our 100 cosmological models can be found in the left panel of Figure 13 (gray curves). We employ the HR suite here and choose to show the HMF at z = 0.55 as a typical redshift of the CMASS galaxies. We also show (red curve) the HMF for the fiducial Planck cosmology (the mean of 28 realizations). The 100 models are taken from the five SLHD slices (hereafter slice 1, 2, ..., 5), each of which consists of 20 different cosmological models, as described in Section 3.1. The variation in the HMF is quite large, and it is not so obvious whether or not the universal form of Equation (32) can explain it.

Zoom In Zoom Out Reset image size

Before constructing an HMF module, we first test the accuracy of the Tinker HMF against our simulation suite (HR runs). Figure 14 compares the simulated HMF with the original Tinker HMF prediction (using the coefficients in Equations (34)–(37)) for each of 20 cosmological models in slice 5 at three redshift bins, z = 0, 0.55, and 1.48. The comparison indicates a larger deviation of the Tinker formula from the simulation results as redshift increases. At z = 0, the ratio of our HR simulation suite to the Tinker MF typically shows 5% scatter from the Tinker MF, with the mean slightly biased from unity (by ∼2%–3%). The overall slope of the HMF is already very well captured by this formula, and the error is mostly on the amplitude. At higher redshifts, both the bias and the scatter grow. At z = 1.48, the slope of the ratio shows a clear mass dependence with a larger bias toward the massive end. We leave further discussion on the inaccuracy of the original Tinker HMF formula to Appendix B, where we discuss that this discrepancy is mainly due to the fact that the Tinker HMF was calibrated against simulations using initial conditions based on the Zel'dovich approximation at an initial redshift that is not high enough.

Zoom In Zoom Out Reset image size

From this exercise, we decide to update some of the parameters in the Tinker HMF. We drop the assumption of the HMF universality and allow the parameters A and a to vary as a function of cosmological models. For the parameters b and c, we keep the values given by Equations (36) and (37), where b determines the slope of the HMF at the low-mass end and c determines the cutoff at the high-mass end. Our HR suite is most accurate over the intermediate range of halo masses, where the overall amplitude and the slope are controlled by the parameters A and a, respectively. Therefore, we recalibrate these parameters for each of our simulations.

The upper middle panel of Figure 13 addresses how the functional form given by Equation (32) can fit the simulated HMF for each of 100 cosmological models with the updated parameters. Here we estimate the best-fitting parameters of a and A that reproduce the simulated HMF. In the fitting, we consider two sources of errors. For high-mass bins, we consider statistical uncertainties due to the Poisson noise of the number counts. We also include a phenomenological penalty term at the low-mass bins, where the correction (Equation (27)) plays a significant role. That is, in the model fitting, we include the following statistical errors in the number counts of halos at each mass bin:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{N}_{{\rm{h}}}}{{N}_{{\rm{h}}}}=\displaystyle \frac{1}{\sqrt{{N}_{{\rm{h}}}}}+\displaystyle \frac{1}{{N}_{{\rm{p}}}},\end{eqnarray} \tag{ 39 }$$

where N_h is the number of halos at each mass bin, and N_p is the number of N-body member particles, M/m_p, at the logarithmic bin center. Moreover, we do not include any mass bin of halos that are defined by N_p < 200. The figure shows that our model HMF generally gives a very good fit to the simulated HMF for each of the 100 cosmological models, where the ratio is very close to unity well within the ±5% accuracy denoted by the horizontal dotted lines. At the high-mass end, the simulation data points are dominated by Poisson noise due to a too-small number of halos per bin. The circle and error bar at each mass bin are the average and scatter in the ratios of the best-fit HMF to the simulated HMF for the 100 cosmological models. For comparison, the red shaded region denotes scatters among the 28 realizations of Planck cosmology, giving an estimate of the sample variance for a volume of 1 (h⁻¹ Gpc)³. The typical accuracy of the model as indicated by the error bars is ∼1% (3%) at 10¹³ (10¹⁴) h⁻¹ M_⊙.

We then compress the data vector, ${\boldsymbol{d}}=({A}_{0},{a}_{0},\,\ldots ,\,{A}_{20},{a}_{20})$, which consists of the fitting parameters A and a at each of 21 redshifts (therefore, 42 data in total) for each cosmological model, using PCA. Combining all data vectors for 100 cosmological models, as well as 28 realizations of the fiducial Planck cosmology (therefore, 128 × 42 = 5376 data points in total), we decompose the data vector for the ith simulation, d_i, into the PCs as

$$\begin{eqnarray}&&{{\boldsymbol{d}}}_{i}=\sum _{j=1}^{n}\,{\alpha }_{i,j}^{\mathrm{HMF}}{{\boldsymbol{e}}}_{j}^{\mathrm{HMF}},\end{eqnarray} \tag{ 40 }$$

where ${{\boldsymbol{e}}}_{j}^{\mathrm{HMF}}$ is the jth eigenvector with 42 components, which is independent of the cosmology or simulation realization, and ${\alpha }_{i,j}^{\mathrm{HMF}}$ is the jth PC coefficient for the ith simulation. After various checks, we find that keeping the six most significant PC coefficients for each cosmological model, corresponding to n = 6 in Equation (40), is sufficient to keep the error induced in this step to a subpercent level. The accuracy level after applying the PCA method is not degraded to an extent easily visible by eye when we compare the lower and upper middle panels of Figure 13. In doing the PCA analysis, we downweight the components A(z) by a factor of 10 compared with a(z) to compensate for their different dynamic ranges.

Our next task is to collect the six PC coefficients (${\alpha }_{i,j}^{\mathrm{HMF}}$ in Equation (40)) for different cosmological models and perform GPR to interpolate each of the six PC coefficients between the sampled cosmological models, each of which is located at a particular position in six-dimensional cosmological parameter space. To do this, we apply the GPR to the PC coefficients for the 80 cosmological models included in slices 1–4 as the training set, excluding the 20 cosmological models in slice 5 (see Section 3.1 for details). Note that we do not include the fiducial Planck cosmology in this GPR either. The upper right panel of Figure 13 compares the GPR HMF with the simulated HMF for each of the 80 cosmological models in the training set. The GP does not perfectly reproduce the results of PCA at each of sampled cosmological models because we take into account the statistical uncertainties of the training data in the regression. Nevertheless, the plot shows that after applying the GPR, the rms in the ratios among different cosmological models are kept below ∼1% (3%) on $M\lesssim {10}^{13}$ (10¹⁴) h⁻¹ M_⊙.

The lower right panel of Figure 13 is the most important plot that gives an assessment of the performance of the HMF emulator for an arbitrary cosmological model. The plot compares the GP-interpolated HMF with the simulated HMF for each of the 20 cosmological models in slice 5 that are not used in the GPR and serve as a cross-validation sample. The HMF emulator achieves great accuracy to predict the HMF, better than a few percent, in the amplitude up to halo masses of a few times 10¹⁴ h⁻¹ M_⊙. The performance is degraded for more massive halos, but the inaccuracy (averaged value denoted by the circle at each bin) is comparable with the statistical scatter. The good performance suggests that our GP method does not suffer from an overfitting to the training set.

While the performance of the emulator can be assessed fairly precisely for halo masses up to ∼10¹⁴ h⁻¹ M_⊙, the scatter among the models appears to be large for cluster-sized halos. While the large scatter could be simply due to the inaccuracy of the emulator, it can be partly due to the large Poisson noise, which can significantly affect the measurements used as the reference due to the small number of available cluster-scale halos in the simulations. As a final check of the accuracy of the emulator at the high-mass end, we compare the emulator prediction to the measurement from the LR simulations, which have bigger volume and are thus less affected by the Poisson noise.

In Figure 15, we show the ratio of the mass function measured from these LR simulations to the emulator prediction, which is trained based on the HR simulation suite. We show the mean and scatter among the 20 cosmological models in slice 5, with individual cosmology result (gray solid curve). We apply the correction of the mass based on the number of member particles (Equation (27)) for the upper set of curves, while we do not apply this to the lower set of curves. Since the size of the correction is pretty large for a mass less than several times 10¹³ h⁻¹ M_⊙ for these sets of simulations, we cannot derive a clear conclusion for these masses given the empirical nature of the correction. For more massive halos, the simulation results after the correction are very close to the emulator prediction. The scatter among the models is much smaller than what we can see in Figure 13, suggesting that the large scatter in the previous figure is indeed mainly due to the large Poisson noise in the reference simulations. However, a closer look at Figure 15 reveals that the mean of the ratio among the cosmological models for cluster-size halos is systematically above unity by ∼1%–2%. This would be a slight inaccuracy of the emulator due to the use of the HR simulation suite with a smaller volume. Since the size of this systematic error is comparable to those discussed for less massive halos based on the HR suite, we do not further consider the possibility of combining the measurements from the LR and HR simulation suites and instead stick to the emulator build based only on the HR suite for the HMF.

Zoom In Zoom Out Reset image size

We now discuss the halo–matter cross-correlation function. We first measure the cross-correlations for 13 mass threshold halo samples with different number densities at each of 21 redshifts from each simulation run, where we define the halo samples in 13 logarithmically spaced bins in the range of ${n}_{{\rm{h}}}=[{10}^{-8.5},{10}^{-2.5}]{({h}^{-1}\mathrm{Mpc})}^{-3}$ (i.e., two bins per decade). For each measurement, we have 140 separation bins (40 logarithmic bins from 0.01 to 5 h⁻¹ Mpc for the direct pair-counting method and 100 linear bins from 5 to 500 h⁻¹ Mpc for the FFT method; see Section 3.6 for details). We thus have 13 × 21 × 140 = 38,220 data points per simulation.

The left panel of Figure 16 shows variations in the halo–matter cross-correlations for 101 different cosmological models in the HR simulation suite for the mass threshold halo sample with a number density of ${10}^{-4}\,{({h}^{-1}\mathrm{Mpc})}^{-3}$ and at z = 0.55. This halo sample is chosen to roughly mimic the number density and large-scale bias of the CMASS galaxies. The plot displays rich cosmological dependences over scales ranging through the one-halo, two-halo terms to BAO scales. We then reduce the dimensionality of the data vector by first resampling the separation bins and then applying a PCA. The former is done using a cubic spline interpolation of the original data points up to 100 h⁻¹ Mpc, with more data points around the one- and two-halo transition scale (i.e., around a Mpc scale), as depicted by the vertical dotted lines in the middle panel of Figure 16. We take 66 data points after the resampling. This procedure degrades the accuracy on small scales by no more than 3% (≲40 h⁻¹ Mpc). One might notice a small wiggly feature around 6 h⁻¹ Mpc. This is due to the grid effect of the FFT method around the switching point to the direct pair-counting method, as described in Figure 3. Our spline function tries to remove this spurious feature to some extent by forcing the curve to be smooth. Since the raw simulation measurements employed here as the numerator still suffer from this artifact, the feature is still present in the ratio.

Zoom In Zoom Out Reset image size

Our data vector still has 13 × 21 × 66 = 18,018 components per simulation, which is quite large. We apply PCA to this data vector. As in the case for HFM, we combine all of the data vectors from 128 simulations (100 for the varied cosmological models plus the 28 Planck cosmology simulations), corresponding to 128 × 18,018 = 2,306,304 data points, and parameterize the halo–matter cross-correlation into its PCs as

$$\begin{eqnarray}&&{\left.{\xi }_{\mathrm{hm}}(x,{n}_{{\rm{h}}},z)\right|}_{i}=\sum _{a=1}^{n}\,{\alpha }_{i,a}^{\mathrm{CCF}}{e}_{a}^{\mathrm{CCF}}(x,{n}_{{\rm{h}}},z),\end{eqnarray} \tag{ 41 }$$

where ${\left.{\xi }_{\mathrm{hm}}(x,{n}_{{\rm{h}}},z)\right|}_{i}$ is the halo–matter cross-correlation at separation x for the halo sample with number density n_h and at redshift z in the ith simulation; ${e}_{a}^{\mathrm{CCF}}$ is the ath PC eigenvectors given as a function of separation x, n_h, and z (cosmology-independent); and ${\alpha }_{i,a}^{\mathrm{CCF}}$ is the ath PC coefficient for the ith simulation. Thus, the eigenvectors $\{{e}_{a}^{\mathrm{CCF}}\}$ describe dependences of the halo–matter cross-correlation on separation, halo sample (halo number density), and redshift. The different eigenvectors, ${e}_{a}^{\mathrm{CCF}}$ and ${e}_{b}^{\mathrm{CCF}}$ with $a\ne b$, are orthogonal to each other. The PC coefficients $\{{\alpha }_{i,a}^{\mathrm{CCF}}\}$ describe the dependences on different simulations, i.e., cosmological models. In applying this PCA, we adopt the weight for each data point by ${n}_{{\rm{h}}}\,x\,[1-\exp (-{x}^{2})]$ (x is in units of h⁻¹ Mpc), such that the data points containing the halo sample of higher number density are upweighted, and the data points at large separation (x) are relatively upweighted; we empirically find the functional form of weight that satisfies these conditions. As shown in the upper right panel of Figure 16, we find that keeping the five significant PC coefficients reproduces the simulation results within a 5% accuracy on x ≲ 30 h⁻¹ Mpc for each of 100 cosmological models, despite the huge dimensionality reduction (from 18,018 to 5). While we can see a relatively large deviation from unity at around 1 h⁻¹ Mpc for a few cosmological models, the rms among the 100 models, as shown by the error bars, is typically below the 2% level on small scales up to ∼20 h⁻¹ Mpc and reaches to 5% at ∼40 h⁻¹ Mpc. Since the halo–matter cross-correlation is a smooth function anyway, the PCA decomposition of the halo–matter cross-correlations works very well.

However, note that the PCA description cannot well describe the correlation at very large separations. The scatter of the curves around unity is consistent with the statistical error of the simulation data due to a finite number of simulation realizations or, equivalently, a finite simulation volume, as shown by the red shaded region, which is estimated from the 28 simulations at the fiducial Planck cosmology. We will instead use a different prescription, the propagator method, for the very large scale to overcome both sample variance in the simulation data and inaccurate modeling, as we will describe later.

We then perform GPR to model the cosmology dependence of the PC coefficients for 80 cosmological models in slices 1–4, which is our training set, similar to Figure 13. Again, note that we do not include the 28 Planck cosmology realizations for this GPR. The upper right panel of Figure 16 shows that the GPR does not degrade the accuracy compared to the results after the PCA by more than 1% for the 80 cosmological models. The lower right panel gives a validation of our emulator; the GP interpolation reproduces equally well the cross-correlation function for each of 20 validation cosmological models in slice 5, which are not used in the GPR. The scatter is similar to the statistical error denoted by the shaded region (the sample variance of 28 Planck cosmology simulations). Thus, again, our GPR interpolation performs well without significant overfitting.

In Appendix F, we show the performance of the GP interpolation for different redshifts, as well as a different sample of halos that are characterized by the different number density in each cosmological model.

We next discuss the halo autocorrelation functions. In this case, we generally need to consider two mass threshold halo samples of different number densities, and thus the dimension of the input data vector is even larger than that for the halo–matter cross-correlation functions. However, we cannot obtain a meaningful signal for halo samples with number densities lower than $\sim {10}^{-6}{({h}^{-1}\mathrm{Mpc})}^{-3}$, unlike the halo–matter cross-correlation function due to the large Poisson noise. Hence, we consider here only eight bins with a high number density out of the 13 bins that we considered for the cross-correlation function. We thus consider 36 (=8(8+1)/2) combinations for the two halo samples to form a matrix of autocorrelation functions at each of 21 redshifts and at each separation bin per simulation.

The left panel of Figure 17 shows variations in the halo autocorrelation functions for 101 cosmological models for the mass threshold halo sample with number density ${n}_{{\rm{h}}}\,={10}^{-4}\,{({h}^{-1}\mathrm{Mpc})}^{-3}$ and at z = 0.55 that are now computed from the LR simulation suite. The plot shows rich cosmological dependences in the halo autocorrelations over the range of separation scales.

Zoom In Zoom Out Reset image size

We then reduce the dimensionality of the data vector by first resampling the separation bins. We originally have 185 bins for each of the autocorrelation functions, 40 from the direct pair-counting method ($0.1\,{h}^{-1}\mathrm{Mpc}\lt x\lt 10\,{h}^{-1}\mathrm{Mpc}$) and 145 for the FFT method (up to $300\,{h}^{-1}\mathrm{Mpc}$). Since the autocorrelation function has a smooth shape at scales below BAO scale, we can significantly reduce the number of sampling points. As in the halo–matter cross-correlation case, we use a cubic spline interpolation to obtain the new data vector, and the vertical dotted lines in the upper middle panel denote the locations of new sampling points (21 points in total). The error level after this procedure is around 3% for the worst cases and 1%–2% in terms of the rms among the models on scales $1\,\lesssim x/[{h}^{-1}\mathrm{Mpc}]\lesssim 40$. The larger error on smaller scales is due to the significant halo exclusion effect.

Even after the above resampling of separation bins, we still have 21 × 21 × 36 = 15,876 data points per simulation. The next task is to reduce the dimensionality of the data vector using PCA. As in the case of the halo–matter cross-correlations, we combine all of the data vectors from 114 simulations (100 cosmological models plus 14 Planck cosmology simulations), corresponding to 114 × 21 × 21× 36 = 1,809,864 data points, and parameterize the halo autocorrelation into its PCs as

$$\begin{eqnarray}&&{\left.{\xi }_{\mathrm{hh}}(x;{n}_{1},{n}_{2},z)\right|}_{i}=\sum _{a=1}^{n}\,{\alpha }_{i,a}^{\mathrm{ACF}}{e}_{a}^{\mathrm{ACF}}(x,{n}_{1},{n}_{2},z),\end{eqnarray} \tag{ 42 }$$

where ${\left.{\xi }_{\mathrm{hh}}(x,{n}_{1},{n}_{2},z)\right|}_{i}$ is the halo correlation function at separation x for the two mass threshold halo samples with number densities ${n}_{{\rm{h}}}={n}_{1}$ and n₂, respectively, and at redshift z in the ith simulation; ${e}_{a}^{\mathrm{ACF}}$ is the ath PC eigenvector given as a function of separation x, n₁, n₂, and z; and ${\alpha }_{i,a}^{\mathrm{ACF}}$ is the ath PC coefficient for the ith simulation. In applying the PCA, we adopt the weight, given as ${n}_{1}{n}_{2}{x}^{2}$, that is given as a function of the two number densities n₁ and n₂ and the separation x. However, since we use the LR suite here, there is a case that we cannot define a sample of halos with the highest number density, e.g., ${n}_{{\rm{h}}}={10}^{-2.5}\,{({h}^{-1}\mathrm{Mpc})}^{-3}$, depending on redshifts and cosmological models. In such cases, we set the weight for PCA to zero. After some experiments, we find that keeping the eight significant PC coefficients well reproduces the simulation results. To be more quantitative, we can maintain the error level of a few percent with a slight degradation toward the large scales ($\gtrsim 10\,{h}^{-1}\mathrm{Mpc}$), as shown in the lower middle panel of Figure 17. Thus, we made a huge dimensionality reduction of the data points (from $1.8\times {10}^{6}$ to 8). The error induced by this GPR is negligibly small except for very large scales, where we stitch to the propagator-based prescription, as we will describe below.

Next, we perform GPR to model the eight PC coefficients for 80 cosmological models: 80 different cosmological models in slices 1–4. The upper right plot of Figure 17 shows that the GP interpolation reproduces the simulation results for each of 80 cosmological models, which is the training set for the GPR, within 1%–5% rms accuracy, depending on the scale. The lower right panel gives a validation of our emulator; the GP interpolation gives 2%–3% rms accuracy on $1\lesssim x/[{h}^{-1}\mathrm{Mpc}]\lesssim 20$ for the 20 cosmological models in slice 5, which is the validation set.

It is worth noting that there are differences in the GPR results between the halo–matter cross- and halo autocorrelations. First, the small-separation data points have large scatters, which is not seen in the cross-correlation function. This is due to the fact that the autocorrelation function is suppressed on these small scales due to the halo exclusion effect. We eventually expect a zero-crossing of the autocorrelation function at the scale where the halo exclusion effect is dominant, and thus the ratio at the zero-crossing scale becomes large and noisy. Second, the accuracy of the emulator prediction is worse than that for the cross-correlation because of the larger noise of autocorrelation measurements due to the Poisson noise originating from the discreteness of halos.

In Appendix F, we show the performance and validation for the halo autocorrelation functions for halo samples of different number densities and redshifts, as well as the cross-correlations between two halo samples with different number densities.

The final piece of our halo modules is the propagator that describes the large-scale correlation functions around BAO scales. Since the halo–matter cross-correlation function involves the propagators of halo and matter, we here study both. We use the LR suite for this purpose, as the damping behavior of the propagator is known to be mainly due to the large-scale bulk motion. Since the estimation of the propagator is done by using the cross-correlation between the halo (or matter) density field and the linear density field used in setting up the initial conditions of the simulations (i.e., Equation (30)), the measured data do not show discreteness noise, unlike what we see in the halo autocorrelation functions. Therefore, we can measure it to a reasonably accurate precision for all 13 number density bins for the mass threshold samples from 10^−2.5 to 10^−8.5 (h⁻¹ Mpc)⁻³ used for the halo–matter cross-correlation functions.

Since the shape of the propagator is roughly a Gaussian with its width being the rms displacement, as shown in Figure 12, we parameterize it as

$$\begin{eqnarray}&&{G}_{a}(k)=\left({g}_{0}+{g}_{2}\,{k}^{2}+{g}_{4}\,{k}^{4}\right)\exp \left(-\displaystyle \frac{{\sigma }_{{\rm{d}},\mathrm{lin}}^{2}{k}^{2}}{2}\right),\end{eqnarray} \tag{ 43 }$$

where g₀ can be interpreted as the linear bias at the large-scale limit and ${\sigma }_{{\rm{d}},\mathrm{lin}}$ is the linear rms displacement for the matter field at the redshift under consideration, which is computed by the linear module. The other coefficients, g₂ and g₄, are fitting parameters, introduced to capture a departure from the simple Gaussian form. In the above, the subscript a on the left-hand side stands for either m (matter) or h (halo), and the factor g₀ is set to unity for the matter propagator.

We show in Figure 18 our modeling detail of this function for halo samples with number density ${10}^{-4}{({h}^{-1}\mathrm{Mpc})}^{-3}$ at z = 0.55. As before, we show variations in the function for the 101 cosmological models in the left panel, with the red shaded region showing the mean and scatter of this function for the fiducial Planck cosmology (the shaded region is also shown in the other panels, though it is heavily overlapped with other symbols and thus difficult to see). We can see that the propagator is always a simple decaying function of wavenumber, with a strong dependence on cosmology in the amplitude and the typical wavenumber at which the curve is decaying.

Zoom In Zoom Out Reset image size

We then fit the data using Equation (43) and show the residual in the upper middle panel. Here and in the other three panels, we normalize the residual by g₀, which is the low-k limit of this function, to see the importance of the residual relative to the overall amplitude of the function. While we see a small wiggly pattern in the residual, the typical amplitude of this pattern is below a few percent level, which is sufficiently small for our purpose.

At this point, we have three fitting parameters per halo sample, and thus 819(=3 × 13 × 21) data points per simulation for the halo propagator. As before, we reduce the dimensionality by applying the PCA. The data points are approximated by

$$\begin{eqnarray}&&{{\boldsymbol{d}}}_{i}^{\mathrm{PRO}}({n}_{h},z)=\sum _{a=1}^{n}\,{\alpha }_{i,a}^{\mathrm{PRO}}{{\boldsymbol{e}}}_{a}^{\mathrm{PRO}}({n}_{h},z),\end{eqnarray} \tag{ 44 }$$

where ${{\boldsymbol{d}}}_{i}^{\mathrm{PRO}}({n}_{h},z)$ is the vector formed with the three fitting parameters for the halo sample with number density ${n}_{{\rm{h}}}$ and at redshift z in the ith simulation, ${{\boldsymbol{e}}}_{a}^{\mathrm{PRO}}$ is the ath PC eigenvector given as a function of n_h and z, and ${\alpha }_{i,a}^{\mathrm{PRO}}$ is the ath PC coefficient for the ith simulation. In applying the PCA analysis, we adopt the weight, simply given as n_h. As before, since we use the LR runs here, there is a case that we cannot define a sample of halos with the highest number density, e.g., ${n}_{{\rm{h}}}={10}^{-2.5}\,{({h}^{-1}\mathrm{Mpc})}^{-3}$, depending on redshifts and cosmological models. In such cases, we set the weight for PCA to zero. After some experiments, we find that keeping the four most significant PC coefficients well reproduces the simulation results, as shown in the lower middle panel of Figure 18. The extra error induced by the PCA is below the 1% level.

The remaining task is the same as before: train a GPR using the 80 cosmological models in slices 1–4 and validate the results using the remaining 20 models in slice 5 (as well as the fiducial Planck cosmology). Although the variance of the residuals among the models shown in the right panels seems to be somewhat larger than that in the middle panels, the prediction of the GP stays within the ±5% band shown by the horizontal dotted lines. More importantly, the accuracy of the GP for the validation set is not degraded compared to that for the training set, implying that there is no problem of overfitting to the training data.

The validation tests for other halo samples, as well as the matter propagator at various redshifts, can be found in Appendix F. In the current implementation, we model the matter propagator by exactly following the same procedure for halos. The only difference is that we have one less free parameter (g₀ should always be unity for matter), and we need only two PCA components to ensure the accuracy.

One application of our emulator is to make predictions for halos in the mass range of galaxies to clusters. Since our modules that compute halo clustering properties directly predict the signals as a function of the cumulative halo number density, which is discretely sampled every 0.5 dex, it might not be so practically useful as it is. To obtain the predictions of halo correlation functions at a given halo mass, we have to interpolate over the sampled number densities and convert it to the halo mass in a given cosmological model. Hence, if one wants to predict the clustering signals as a function of halo mass, an inaccuracy in the conversion from the mass to the number density can be a new source of error.

In Figure 19, we study how the conversion from the cumulative halo number density to a target halo mass using the HMF module causes a possible error in the predictions of halo correlation functions. To do this, we consider the emulator outputs of halo correlation functions for two mass threshold samples with ${M}_{\min }={10}^{13}$ and ${10}^{14}{h}^{-1}{M}_{\odot }$ for the fiducial Planck cosmology at z = 0.55. Then we use the following method to propagate a possible error in the halo number density into an error in predicting the halo correlation functions as a function of the halo mass. (i) We first use the HMF module to compute the cumulative number density for the halo mass thresholds, ${M}_{\min }={10}^{13}$ and ${10}^{14}{h}^{-1}{M}_{\odot }$. (ii) We multiply the number density by a factor of 0.96, 0.98, 1.02, or 1.04, which is intended to mimic a possible error in the number density calibration by −4%, −2%, 2%, or 4%, respectively. (iii) We then obtain the emulator predictions of halo correlation functions by inserting the shifted values of halo number density in the emulator. Here a ±4% error in the cumulative halo number density is considered as a rather pessimistic case because Figure 13 shows that a typical error in the mass function is smaller in terms of the rms among the models (∼1% (3%) at ${10}^{13}({10}^{14}){h}^{-1}{M}_{\odot }$). In addition, the error in the HMF seen in Figure 13 would be partly canceled when we consider the error on the cumulative halo number density. The upper panel shows the ratios of the shifted halo–matter cross-correlation functions, ${\xi }_{\mathrm{hm}}(x)$, to the fiducial prediction. The figure shows a constant shift in the two-halo regime, a slightly larger shift in the one-halo term, and a bump-like feature at transition scales between the two regimes. These are caused by changes in the linear bias and mass profile, respectively. The size of the fractional shift in the cross-correlation function is smaller than that on the cumulative mass function, with a slight decreasing trend toward higher masses. In the lower panel, the autocorrelation function, ${\xi }_{\mathrm{hh}}(x)$, shows a larger shift in the two-halo regime, reflecting the fact that it scales as bias squared. A sharp feature can be found where the halo exclusion effect kicks in. Since the latter part is dominated by the one-halo term in the case of the galaxy correlation function, the final shift would be much smaller. Even with the pessimistic case of a 4% error in the cumulative mass function, the induced shift in the correlation functions are well within the ±5% band and mostly within the ±3% level. When we consider a realistic error on the cumulative mass function (i.e., a few percent or below), the error on the correlation function arising from this is smaller than the typical error in the emulator in both cases.

Zoom In Zoom Out Reset image size

As another example of the applications, we show in Figure 20 the output of the emulator for the large-scale bias as a function of the halo number density (symbols). Here the large-scale bias is defined as the fitted parameter g₀ in Equation (43), which is the $k\to 0$ limit of the propagator (Equation (28)). The plot shows the result for the fiducial Planck cosmology and at z = 0.55. We interpolate these data points using the cubic spline function to make a prediction at any halo number density in the range shown here. Note that we use the logarithm of the halo number density, instead of the raw values of the number density, for which our sampling is uniform.

Zoom In Zoom Out Reset image size

Once the spline interpolator is ready, we can compute the halo bias as a function of the halo mass or peak height, if one prefers, by first using the HMF module to convert the number density to the minimum halo mass and then taking a finite difference derivative to have the bias at a specific desired halo mass scale. For this derivative, we employ ±1% changes in the mass as the default step size. The dependence of the results on the step size is much weaker than the typical accuracy of the emulator, as shown in Figure 21. The result is shown in Figure 22 as a function of the peak height $\nu \equiv {\delta }_{{\rm{c}}}/{\sigma }_{M}$, with ${\delta }_{{\rm{c}}}=1.686$. Now we show the results for different cosmologies at three different redshifts. We vary ${{\rm{\Omega }}}_{{\rm{m}}}$, keeping the flatness and fixing the present-day amplitude of the linear matter power spectrum, ${\sigma }_{8}$, to its fiducial value. We compare the results with the fitting formula by Tinker et al. (2010), as denoted by the thin solid line, which is independent of redshift or cosmology. The fractional differences from Tinker et al. (2010) at the three redshifts are plotted together in the bottom panel. Our emulator prediction is overall consistent with the fitting formula, with the accuracy no worse than 10% over all the ranges examined here. Such an inaccuracy of Tinker et al. (2010) is also pointed out in the previous work (e.g., Li et al. 2016). We confirm that the bias function is rather universal, with little dependence on cosmology or redshift.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another interesting feature of the bias is its scale dependence around the BAO scale. We implement the same spline interpolation for the fitting parameters g₂ and g₄ in the propagator. This allows us to estimate the propagator at any halo mass. The tree-level calculation, Equation (29), gives us a prediction of the correlation functions, and we already show that the BAO scale is well described by this simple model, as illustrated in Figures 9–11. We show in Figure 23 the square root of the ratio of the halo and matter correlation functions at z = 0. We normalize it by the linear bias factor g₀ such that the ratio becomes unity when the bias is independent of scale. We consider four halo samples with different number densities as written in the figure legend. The result indicates that the BAO peak structure can be boosted for low number density samples (i.e., when only massive halos are included in the sample). This is fully consistent with the expectation by the peak model (compare our results with Figures 7 and 8 of Desjacques et al. 2010). Also, this feature was previously found in numerical simulations (e.g., Angulo et al. 2014; Crocce et al. 2015). This kind of prediction is possible because our model has the freedom to control the damping of the BAO feature in terms of the two free parameters, g₂ and g₄, in the propagator, Equation (43), in addition to the damping due to the typical random displacements of matter, ${\sigma }_{{\rm{d}},\mathrm{lin}}$.

Zoom In Zoom Out Reset image size

We implement the same cubic spline interpolation for other quantities, such as ${\xi }_{\mathrm{hm}}(r;{n}_{h})$ or ${\xi }_{\mathrm{hh}}(r;{n}_{1},{n}_{2})$, with the latter using the bivariate cubic spline for two number densities, n₁ and n₂. We do not find any sizable error originating from this interpolation, as all of the quantities vary rather smoothly with (the logarithm of) the halo number density. Analogously, the redshift dependence is interpolated with the cubic spline function. As is clear from Figure 22, the dependence of bias on redshift is weak, and thus the same interpolation scheme works fine. The situation is similar for the other interpolated quantities. We show in Figures 24 and 25 our interpolation of the halo–matter cross- or halo autocorrelation function over the number density and redshift. We first construct a data matrix at the locations, as depicted by the dots, based on the GP and PCA methods, and then we perform a cubic spline interpolation for each dimension.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the halo–matter cross-correlation function and the propagator are very well determined down to a quite low halo number density, ${10}^{-8.5}{({h}^{-1}\mathrm{Mpc})}^{-3}$—corresponding to very massive halos, thanks to the fact that they are given by the cross-correlation with the matter field—the halo autocorrelation function instead suffers from severe Poisson noise, especially at such a high-mass end. We thus switch to a simple scaling, ${\xi }_{\mathrm{hh}}(r;{n}_{1},{n}_{2})=[{g}_{0}({n}_{1})/{g}_{0}({n}_{\min })]\,{\xi }_{\mathrm{hh}}(r;{n}_{\min },{n}_{2})$, when the number density n₁ is below the minimum number density ${n}_{\min }={10}^{-5.75}{({h}^{-1}\mathrm{Mpc})}^{-3}$. In the above, the bias factor, ${g}_{0}({n}_{i})$, is computed again in the module that computes the propagator (i.e., the function plotted in Figure 20). We do the same when n₂ is below the threshold; we simply multiply the ratio of the large-scale bias one more time. While this might be a reasonable approximation on large scales, it cannot properly reproduce the correlation functions around scales where the halo exclusion effect is not negligible. Nevertheless, our current implementation does not lead to a severe error for a sample of galaxies such as the CMASS sample, because the small-scale correlation function is mainly described by the one-halo term.

Now we have predictions of ${\xi }_{\mathrm{hm}}$, ${\xi }_{\mathrm{hh}}$, and the propagator given as a function of halo number density and redshift within the ranges relevant for the resolution and output redshifts of our simulations. We combine all of these predictions to obtain a well-behaved prediction over a wide range of separations. We do this by smoothly stitching the two predictions as

$$\begin{eqnarray}&&{\xi }_{\mathrm{ab},\mathrm{full}}(x)=D(x){\xi }_{\mathrm{ab},\mathrm{direct}}(x)+\left[1-D(x)\right]{\xi }_{\mathrm{ab},\mathrm{tree}}(x),\end{eqnarray} \tag{ 45 }$$

where "ab" is either "hm" or "hh," and we use the damping function defined as

$$\begin{eqnarray}&&D(x)=\exp \left[-{\left(\displaystyle \frac{x}{{x}_{\mathrm{switch}}}\right)}^{4}\right].\end{eqnarray} \tag{ 46 }$$

We find that ${x}_{\mathrm{switch}}=60\,{h}^{-1}\mathrm{Mpc}$ provides a reasonably good model for both the auto- and cross-correlation functions over the range of scales demonstrated in Figure 26.

Zoom In Zoom Out Reset image size

In summary, the current implementation of our emulator works in a parameter sampler as follows. When a new cosmological model is proposed, the code first calls a GP interpolator to evaluate the coefficients for the PCs for all of the statistics we consider here. Then, combining these coefficients with the eigenvectors, it computes the statistics at all of the redshifts and mass and separation bins to form a data table. One function call of our high-level interface to set the cosmological parameters does all the tasks up to here internally. Now, a user can further call other high-level functions prepared for each statistics. These functions accept a redshift, a halo mass (either a threshold mass or a target mass scale), and a set of separations at which the correlation function should be evaluated. In this final step, the code finds the values by calling a spline interpolator over the table created in the previous step. For users who wish to compute galaxy statistics, a separate module can be used with additional parameters describing the HOD model. This module internally calls the functions for the halo statistics and integrates them over the halo mass with the product of the mean HOD and the HMF as a weight. We also prepare functions computing projected statistics, which works similarly.

Now we can compute the three main halo statistics—the HMF, halo–matter cross-correlation function, and halo autocorrelation function—for an arbitrary cosmological model that is covered by our sampled cosmological models within the flat wCDM cosmologies. Using the results, we can obtain how these halo statistical quantities vary with cosmological parameters, as demonstrated in Figure 1, which gives their dependences on ${{\rm{\Omega }}}_{{\rm{m}}}$.

We can further predict in detail, for instance, the density profile of dark matter halos from our emulator. Properties and cosmological dependences of the mass density profiles around halos have been extensively studied (e.g., Navarro et al. 1996). The emulator output of the halo–matter cross-correlation on small separations can be used to study the mass density profiles for halos whose masses are in the range supported by the emulator.¹⁷

In Figure 27, we compare the profiles from the emulator (solid) with the best-fit Navarro–Frenk–White (NFW) profiles (dashed) for halos of different masses. For the fitting, we included the data over the range of radii from twice the softening scale to 80% of R_200m. The middle panel shows that the NFW profile gives a good fit to a fractional accuracy better than about 5% up to the virial radius (R_200m), beyond which the NFW profile no longer reproduces the simulation results.

Zoom In Zoom Out Reset image size

Using the NFW fitting results, we can also study how the concentration parameter varies with halo mass, as well as cosmological models over different redshifts, where we define the concentration parameter, ${c}_{200{\rm{c}}}$, by the ratio of the radius within which the density is 200 times the critical density to the scale radius determined by the NFW fit. Figure 28 shows ${c}_{200{\rm{c}}}$ as a function of the peak height, $\nu ={\delta }_{c}/{\sigma }_{M}$, for cosmologies with different ${{\rm{\Omega }}}_{{\rm{m}}}$ at three redshifts. Similar to the previous plots, we keep the spatial flatness and vary the normalization A_s such that σ₈ is kept unchanged for models with different Ω_m. While the relation seems to be universal at high redshift with little dependence on Ω_m, we can see a clear dependence at lower redshifts. The increasing trend of c_200c as decreasing redshift, as well as its positive Ω_m dependence, can be found in Diemer & Kravtsov (2015). Further study of the dependence of the concentration–mass relation on the cosmological parameters can be found in Kwan et al. (2013).¹⁸ In this way, our emulator approach automatically incorporates a possible nonuniversality of the concentration–mass relation. This is quite different in the standard analytical halo model approach, where one usually employs a simulation-calibrated scaling relation for the concentration. We would also like to note that such a calibration of the concentration is often done for a specific cosmological model.

Zoom In Zoom Out Reset image size

Now we focus on the halo mass density profiles at radii larger than ${R}_{200{\rm{m}}}$ in Figure 27, where NFW no longer gives a good fit. The figure shows a clear feature of the transition from the one-halo to the two-halo regime. This feature recently drew attention as a possible "physical" outer boundary of a halo associated with the first orbital apocenter of accreted matter after its infalling, dubbed the "splashback" feature (Diemer & Kravtsov 2014; also see More et al. 2016, for the first detection from observational data). This feature has already been studied with our Dark Quest simulation suite in Okumura et al. (2018a, 2018b), with particular attention to the feature in the velocity statistics around halos.

The feature can be found from the bottom panel of Figure 27, where we show the logarithmic slope of the mass density profile. We obtain this by first fitting the emulator results by the functional form proposed in Diemer & Kravtsov (2014; hereafter the DK fit), and then we take the derivative. We do this for the separation range again from twice the softening scale but to four times the radius R_200m to cover both the one- and two-halo regimes. The best-fit model is shown in the top panel by the dotted lines (but they are difficult to distinguish from the solid lines; they are almost on top of each other), and the ratio of the emulator results is shown in the middle panel. The accuracy of the fit is similar or better than the NFW form, and it remains to be a good fit to much larger scales. Now we can see in the bottom panel that the derivative based on the DK fit shows a sharp dip with a slope steeper than the outer NFW slope (i.e., −3), marking the location of the splashback radius.

Figure 29 shows the splashback radius, R_sp, for various cosmological models and redshifts, where we define R_sp by the location of the minimum logarithmic slope of the DK fit. For clarity, we plot the ratio, ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$, as a function of the peak height. Overall, R_sp is similar to ${R}_{200{\rm{m}}}$, with a slight decreasing trend as a function of the peak height. In addition, the ratio is higher for cosmological models with larger Ω_m. These trends are in qualitative agreement with the fitting formulae in More et al. (2015a), in which the dependence is encoded in the redshift-dependent density parameter ${{\rm{\Omega }}}_{{\rm{m}}}(z)$ in addition to the peak height ν (see also Adhikari et al. 2014)¹⁹

Zoom In Zoom Out Reset image size