The Mira-Titan Universe. III. Emulation of the Halo Mass Function

Matter in the universe is known to cluster in the form of localized, clumpy distributions called halos. The abundance of massive halos as a function of total mass (dark matter and baryons) and redshift—the mass function—depends sensitively on cosmological parameters. In the context of this paper, the term "massive" refers to halo masses characteristic of galaxy clusters, a rich and well-studied class of objects, with a central place in modern cosmology (Mulchaey et al. 2004). Theoretical predictions of the mass function on observationally relevant mass scales—an essentially nonlinear quantity—can now be carried out with good accuracy using large-scale cosmological simulations. Measurements of cluster-scale halo masses can then be used to probe the evolution of cosmic structure growth and to constrain cosmological parameters (Holder et al. 2001; for a review, see Allen et al. 2011). In practice, such halos are detected via the (localized) presence of galaxies (Gladders & Yee 2000; Rykoff et al. 2014), gravitational lensing signatures (Miyazaki et al. 2018), and/or the presence of hot ionized gas contained within the deep gravitational potential wells characteristic of rich galaxy groups and galaxy clusters (via X-ray emission (Böhringer et al. 2001; Burenin et al. 2007; Pacaud et al. 2016) and Sunyaev–Zel'dovich effect observations (Bleem et al. 2015; Planck Collaboration et al. 2016; Hilton et al. 2018)). The cosmological constraints resulting from the analyses of cluster samples already allow for competitive measurements of the growth history of the universe (Mantz et al. 2015; Bocquet et al. 2019; Costanzi et al. 2019; Zubeldia & Challinor 2019). In the near future, large samples of thousands of clusters, combined with a mass calibration that is accurate at the few-percent level will allow for greatly improved measurements of the cluster mass function (e.g., Abbott et al. 2020).

From the theoretical perspective, halos suffer from the lack of a strict definition, and therefore, the notion of "halo mass" inherits a certain ambiguity (see, e.g., White 2001). Aside from this fact, the halo mass itself is not a direct observable. For these reasons alone, connecting observations to a theoretically obtained mass function is a nontrivial task. In the absence of a robust forward-modeling approach based on detailed simulations, the current state of the art relies on empirical "mass–observable" relations to connect theory and observations. In such an approach, one important limitation is the accuracy with which the theoretical mass function is known (for a given definition of mass), as the relevant cosmological parameters are varied. This general topic is the focus of the work presented here.

Early estimates of the mass function applied the spherical collapse model to the linear matter density field (Press & Schechter 1974), further formalized via the excursion set approach (Bond et al. 1991). In this approximate methodology, all dependence of the mass function on redshift and cosmology is completely described by the rms fluctuations σ(M, z) in the linear matter power spectrum P(k, z). It is not obvious, however, that this "universal" prediction of the mass function would be sufficiently accurate. Initial results from numerical studies of the mass function based on N-body simulations found that this universality did hold at an approximate level, given a certain halo mass definition (Jenkins et al. 2001). A number of fits for the mass function have since been derived from N-body simulations, based on the idea that the universal form can be applied to multiple cosmological models across a wide range of redshifts (see, e.g., Sheth & Tormen 1999; Jenkins et al. 2001; Springel et al. 2005; Heitmann et al. 2006b; Warren et al. 2006). (For a discussion on the possible halo definition dependence of the universal form, see White 2002.)

Despite the remarkable success of the approach described above, as simulation results were further refined, it was discovered that the redshift evolution of the mass function, even for ΛCDM, deviates from the universal prediction (at the 5%–10% level) and that this deviation had to be explicitly fitted (e.g., Lukić et al. 2007; Reed et al. 2007; Cohn & White 2008; Tinker et al. 2008; Crocce et al. 2010; Bhattacharya et al. 2011; Courtin et al. 2011). Furthermore, it was found that the universal mass function fit established for a reference ΛCDM model can only be extrapolated to wCDM models (with w ≃ −1) at about 10% accuracy (e.g., Bhattacharya et al. 2011).

In order to obtain precision cosmological constraints from cluster samples, one therefore needs to proceed beyond the use of fitting functions of the type discussed above, given their limitations in accuracy and parametric coverage. Reducing the systematic uncertainty in the theoretical prediction by roughly an order of magnitude is not an easy task (for a detailed discussion of these issues, see, e.g., Bhattacharya et al. 2011), especially if one wishes to include the effects due to nonzero neutrino mass and dynamical dark energy, and also account for the influence of baryons, all of which have potential ramifications for the behavior of the mass function.

In the cosmological context, a direct numerical approach to this problem with a finite number of sufficiently accurate simulations is in fact possible using a combination of efficient sampling strategies, Bayesian statistical methods, and machine-learning-based data reduction and interpolation—a process termed emulation (e.g., Heitmann et al. 2006a; Habib et al. 2007; Higdon et al. 2010). The end result of the emulation process, an emulator, is an oracle that, given a set of input parameters, yields an essentially instantaneous prediction for a set of summary statistics, with well-defined errors.

To construct an emulator, numerical simulations are first run for a set of cosmologies sampling a bounded parameter space. Statistical data reduction and machine-learning-based interpolation techniques then yield the desired predictions for observables for any set of parameters contained within the sampled region. Emulators have been used to successfully predict the mass function and other nonlinear summary statistics characteristic of cosmological structure formation such as the nonlinear matter power spectrum, halo bias, and halo concentration (Lawrence et al. 2010, 2017; Kwan et al. 2013, 2015; Heitmann et al. 2014; McClintock et al. 2019; Nishimichi et al. 2019).

In this paper, we present a mass function emulator using the Mira-Titan Universe suite of N-body simulations (Heitmann et al. 2016). This simulation suite includes the effects of massive neutrinos as well as dynamical dark energy and is therefore suited for current and next-generation cosmological surveys that are searching for deviations from ΛCDM. The box sizes (2.1 Gpc) and mass resolution (3200³ particles) of the simulations were designed to measure the mass function at mass scales starting at $\sim {10}^{13}\,{M}_{\odot }/h$ with good resolution of individual halos (at least 1000 particles per halo) as well as good statistics of halo masses in individual mass bins on group/cluster mass scales. Our emulator provides the mass function up to redshift z = 2, with percent-level accuracy at group-scale masses and better than ∼10% accuracy at ${10}^{15}{M}_{\odot }/h$. This new emulator is a complement to the Mira-Titan emulator for the matter power spectrum presented in Lawrence et al. (2017).

This paper is structured as follows. The Mira-Titan Universe simulations and the extraction of binned halo catalogs are discussed in Section 2. We describe the construction of the halo mass function emulator and verify its accuracy in Section 3. We discuss existing predictions for the mass function in Section 4. In Section 5, we use the Mira-Titan Universe simulations to discuss the limits of the cosmological universality of the mass function. We conclude with a summary and outlook in Section 6.

The simulations used in this work were run using the Hardware/Hybrid Accelerated Cosmology Code (HACC) cosmological N-body code (Habib et al. 2016). The suite of simulations is named the Mira-Titan Universe after the supercomputers used to produce it (the IBM BlueGene/Q system Mira at Argonne and the GPU-accelerated system Titan at Oak Ridge) and is introduced in Heitmann et al. (2016). Here, we summarize the relevant aspects and refer the reader to the original paper and references therein for further details.

2.1. Cosmological Modeling: $\nu {w}_{0}{w}_{a}$ CDM

The Mira-Titan Universe is a suite of cosmological N-body simulations that are realizations of 111 different cosmologies, which we refer to as M001–M111. The cosmological parameters are drawn from an eight-dimensional parameter space {${{\rm{\Omega }}}_{{\rm{m}}}{h}^{2}$, ${{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}$, ${{\rm{\Omega }}}_{\nu }{h}^{2}$, ${\sigma }_{8}$, h, n_s, ${w}_{0}$, ${w}_{a}$}. All models are spatially flat (Ω_k = 0). We consider a dynamical dark energy equation of state, using the common parameterization (Chevallier & Polarski 2001; Linder 2003)

$$\begin{eqnarray}&&w(a)={w}_{0}+{w}_{a}(1-a).\end{eqnarray} \tag{ 1 }$$

We further define the parameter combination

$$\begin{eqnarray}&&{w}_{b}\equiv {\left(-{w}_{0}-{w}_{a}\right)}^{1/4}\end{eqnarray} \tag{ 2 }$$

as we will not consider the (w₀, w_a) parameter space but (w₀, w_b) (see discussion after Equation (10)).

Massive neutrinos are treated at the background level instead of being simulated as a separate particle species, essentially an expansion in the neutrino mass fraction, f_ν, in the spirit of Saito et al. (2009). A detailed description of our approach is given in Upadhye et al. (2014) and Heitmann et al. (2016). Its validity with regard to power spectrum measurements on large to quasilinear scales is discussed extensively in Upadhye et al. (2014) and the effect of neutrinos on the mass function is investigated for two models in Biswas et al. (2019), including one for high-mass neutrinos. For completeness, we provide a short summary of the approach here. As mentioned above, we do not include the nonlinear evolution of the neutrinos explicitly in the simulation but rather evolve the cold dark-matter-baryon component only. The neutrinos are included in the background evolution (and the initial condition) and therefore do affect matter clustering in the nonlinear regime. However, neutrino clustering itself is not taken into account. Particular care is given to the setup of the initial conditions. We include the neutrino contribution in the transfer function and generate a linear power spectrum at z = 0 with the chosen σ₈. We then move the linear power spectrum back to the initial redshift using the scale-independent growth function (note that the growth function including the nonlinear neutrino evolution is scale dependent). The growth function includes all species in the homogeneous background. The use of the scale-independent growth function is important because the evolution is carried out only for the dark-matter-baryon component. This approach ensures that the fully developed total (clustered) matter power spectrum at redshift z = 0 is correctly normalized to a given σ₈. In this approximate approach, the linear fluctuations in the neutrinos have to be added separately to get the power spectrum, where the nonlinear contribution of the neutrinos is ignored. At the current observationally favored neutrino mass range, ∼0.1 eV or less, the nonlinear clustering of neutrinos has a negligible effect on halo masses, while at the very upper end of the masses considered here (∼1 eV), the effect on cluster-scale halos is subdominant compared to the overall suppression of the mass function due to neutrino free streaming and is of the order of the overall accuracy of the emulator.

The simulation suite used in this work covers gravity-only simulations. The expectation is that direct modeling of baryonic effects (as well as approaches based on postprocessing gravity-only simulations) will be added over time (Heitmann et al. 2016).

2.2. The Design

The Mira-Titan Universe was specifically designed to produce simulation data for emulators. The volume of the sampled parameter space is therefore a compromise between covering wide parameter ranges and the sparsity this implies. The sampling design underlying the Mira-Titan Universe has an in-built notion of sequential convergence and is aimed at generating multiple emulators at percent-level error levels, including the halo mass function (as demonstrated in a first test in Heitmann et al. 2016).

The cosmological parameters are chosen within the following boundaries:

$$\begin{eqnarray}&&0.12\leqslant {{\rm{\Omega }}}_{{\rm{m}}}{h}^{2}\leqslant 0.155\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&0.0215\leqslant {{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}\leqslant 0.0235\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&0.0\leqslant {{\rm{\Omega }}}_{\nu }{h}^{2}\leqslant 0.01\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&0.7\leqslant {\sigma }_{8}\leqslant 0.9\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&0.55\leqslant h\leqslant 0.85\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&0.85\leqslant {n}_{s}\leqslant 1.05\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&-1.3\leqslant {w}_{0}\leqslant -0.7\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&0.3\leqslant {w}_{b}\leqslant 1.3.\end{eqnarray} \tag{ 10 }$$

The cosmological hypercube is spanned by w_b instead of w_a to ensure a better coverage of models with ${w}_{0}+{w}_{a}$ close to zero (Heitmann et al. 2016). As a result, w_a is jointly constrained with w₀; the smallest allowed value is w_a = −2.16 and the highest allowed value is w_a = 1.29. The range in ${{\rm{\Omega }}}_{\nu }{h}^{2}$ corresponds to a range in the sum of neutrino masses $0\leqslant \sum {m}_{\nu }\leqslant 0.94\,\mathrm{eV}$. While the upper limit is significantly higher than current results for a ${\rm{\Lambda }}\mathrm{CDM}$ background cosmology, the constraint on ${{\rm{\Omega }}}_{\nu }{h}^{2}$ significantly degrades when w₀ and w_a are also allowed to vary (see Figure 1).

Zoom In Zoom Out Reset image size

The choice of the 111 design models and their space-filling properties are discussed in Section 3 in Heitmann et al. (2016). Note that a tessellation-based nested design strategy was specially developed in order to obtain sequential convergence, allowing reuse of previously run simulations: the first 26 models (M011–M036) represent a complete design in their own right, which is further refined when adding the next 29 models to obtain a design with 55 models. Adding yet another 46 models gives a full design with 101 models. To improve coverage of the edge of parameter space where $\sum {m}_{\nu }=0$, another set of 10 models (M001–M010) with massless neutrinos are chosen according to a seven-parameter symmetric Latin hypercube design. All cosmologies in the Mira-Titan design are listed in Table A1 and are visualized in Figure 1.⁷ As discussed earlier, the design hypercube is constructed using w₀ and w_b, which leads to the impression of uneven sampling in (w₀, w_a) space.

For each of the 111 design cosmologies, a 2.1 Gpc box simulation was run evolving 3200³ particles. The particle masses vary between 7.23 × 10⁹ M_⊙ and 1.22 × 10¹⁰ M_⊙ depending on the simulated cosmology, comfortably resolving group and cluster-scale halos. Key features of the simulations are summarized in Table 1.

Table 1.  Mira-Titan Universe Simulations Used to Construct the Emulator

Model	Box Size	Nparticle	Force Res.	Mparticle	$\min ({M}_{200{\rm{c}}})$
M001–M111	2.100 Gpc	32003	6.6 kpc	$1.08\times {10}^{10}\left(\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{m}}}-{{\rm{\Omega }}}_{\nu }}{0.28}\right){\left(\tfrac{h}{0.7}\right)}^{2}{M}_{\odot }$	${10}^{13}{M}_{\odot }/h$
M006	2.091 Gpc	32003	6.6 kpc	1.02 × 1010M⊙	${10}^{13}{M}_{\odot }/h$
M023	2.085 Gpc	32003	6.6 kpc	9.96 × 109M⊙	${10}^{13}{M}_{\odot }/h$
M046	1.865 Gpc	32003	6.6 kpc	7.23 × 109M⊙	${10}^{13}{M}_{\odot }/h$

Note. The mass cuts $\min ({M}_{200{\rm{c}}})$ ensures we only use well-resolved halos with N_particle > 1000. The simulations for models M006, M023, and M046 were accidentally run with slightly smaller box sizes than 2.1 Gpc. As we build an emulator for the spatial halo number density, it is trivial to account for this change in volume.

Download table as:  ASCIITypeset image

An important verification test of the emulator accuracy is to compare its mass function prediction directly with numerical results from additional (off-design) simulations that were not used for the construction of the emulator. For this purpose, a set of four additional 2.1 Gpc box simulations was run (T001–T004). Finally, another 2.1 Gpc simulation box was run using a fiducial ΛCDM cosmology (M000) that has also been used for other, larger simulations (e.g., the Outer Rim simulation, Heitmann et al. 2019). All of the five reference simulations evolve 3200³ particles as in our main simulation suite.

2.3. Halo Catalogs and Binned Halo Mass Function

We first identify halos using a fast parallel friends-of-friends (FOF) halo finder with link length set to b = 0.168. Spherical overdensity (SO) halo catalogs are then built out from the potential minimum of each FOF halo. (In this second process, we use all of the local simulation particles, not just those found by the FOF algorithm.) SO halo masses are specified via the overdensity criterion Δ = 200ρ_crit. We create halo catalogs for eight snapshots that are roughly equally spaced in time within 0 ≤ z ≤ 2.02. The redshifts are z ∈ {0.00, 0.10, 0.24, 0.43, 0.66, 1.01, 1.61, 2.02}. Throughout the simulation suite, we only consider halos above ${10}^{13}{M}_{\odot }/h;$ this cut ensures that all halos are well resolved with at least 1000 particles.

For each halo catalog, we apply a binning in mass which is a compromise between good mass resolution at low mass where halos are most abundant and a sufficiently large number of halos per bin at high mass. First, we bin all halo catalogs with a default bin width ${\rm{\Delta }}{\mathrm{log}}_{10}M=0.1$. Then, we combine the highest-mass bins such that the resulting (wide) high-mass bin contains at least 20 halos. Using 4³ spatial jackknife subvolumes, we compute covariance matrices between mass bins to account for noise in the binned halo catalogs due to sample variance and shot noise. We checked that our mass function emulator is robust to changes in the details of the binning scheme and the estimation of the statistical noise. We ignore correlations between mass bins across different time snapshots. In Figure 2, we show the binned halo mass functions for the 111 design cosmologies.

Zoom In Zoom Out Reset image size

In this section, we discuss our emulator framework. We treat each snapshot separately and thus construct eight independent emulators. The key steps in building the emulators are summarized below; we discuss each of these steps in detail in the following subsections.

• 1.  
  For every input cosmology, obtain a smooth halo mass function model from the halo catalog.
• 2.  
  Perform a principal component analysis (PCA) on the (logarithm of the) smooth mass function models.
• 3.  
  Keep the first four eigenvectors as basis functions.
• 4.  
  Fit the four basis functions to each halo catalog to obtain a four-dimensional posterior parameter distribution for each model.
• 5.  
  Set up Gaussian process (GP) regression to interpolate between the four fit parameters, accounting for their covariances.
• 6.  
  Set the hyperparameters of the GP.
• 7.  
  Verify the emulator predictions using hold-out tests and spot checks against additional simulations.

The emulator is now complete and, for a given requested cosmology, returns the halo mass function and an error estimate.

For our emulator and its construction, we consider mass in units of ${M}_{\odot }/h$ and volume in units of ${\left(\mathrm{Mpc}/h\right)}^{3}$. With this convention, the conversion from unit volume to the volume contained in a survey solid angle does not explicitly depend on the Hubble constant h.

3.1. Constrained Piecewise Polynomial Halo Mass Function Fits

To ensure that the final emulator predictions for the mass function are smooth, we need to convert the binned halo catalogs into smooth functions of mass. Due to the large variance of the mass function at high mass, we found smoothing kernels not to be an adequate approach. Inspired by the approximate power-law behavior of the mass function at low masses, we instead fit a model consisting of constrained piecewise second-order polynomials in log mass. In each segment i, we have

$$\begin{eqnarray}&&\mathrm{ln}{\left(\displaystyle \frac{{dn}}{d\mathrm{ln}M}\right)}_{i}={a}_{i}+{b}_{i}\mathrm{ln}\left(\displaystyle \frac{M}{{M}_{\mathrm{piv}}}\right)+{c}_{i}\mathrm{ln}{\left(\displaystyle \frac{M}{{M}_{\mathrm{piv}}}\right)}^{2}.\end{eqnarray} \tag{ 11 }$$

We constrain the individual segments to connect smoothly by requiring the function and its derivative to be continuous. For every two consecutive segments joining at mass M_join, this leads to a constraint on the function value

$$\begin{eqnarray}&&\begin{array}{l}{a}_{i}+{b}_{i}\mathrm{ln}\left(\displaystyle \frac{{M}_{\mathrm{join}}}{{M}_{\mathrm{piv}}}\right)+{c}_{i}\mathrm{ln}{\left(\displaystyle \frac{{M}_{\mathrm{join}}}{{M}_{\mathrm{piv}}}\right)}^{2}\\ =\,{a}_{i+1}+{b}_{i+1}\mathrm{ln}\left(\displaystyle \frac{{M}_{\mathrm{join}}}{{M}_{\mathrm{piv}}}\right)+{c}_{i+1}\mathrm{ln}{\left(\displaystyle \frac{{M}_{\mathrm{join}}}{{M}_{\mathrm{piv}}}\right)}^{2}\end{array}\end{eqnarray} \tag{ 12 }$$

and its derivative

$$\begin{eqnarray}&&{b}_{i}+2{c}_{i}\mathrm{ln}\left(\displaystyle \frac{{M}_{\mathrm{join}}}{{M}_{\mathrm{piv}}}\right)={b}_{i+1}+2{c}_{i+1}\mathrm{ln}\left(\displaystyle \frac{{M}_{\mathrm{join}}}{{M}_{\mathrm{piv}}}\right).\end{eqnarray} \tag{ 13 }$$

In addition, we require

$$\begin{eqnarray}&&{c}_{i}\lt 0\end{eqnarray} \tag{ 14 }$$

to encapsulate our intuition that the mass function becomes steeper for increasing mass.

We define segments with constant length in log mass with four segments per decade in mass. We define the fourth segment as the pivot segment for which we define the overall amplitude parameter a₄ as well as b₄ and c₄. With these parameters set, the only free parameters left are the c_i for all other segments; the remaining a_i and b_i are set according to the constraints from Equations (12) and (13). We define ${M}_{\mathrm{piv}}={10}^{14}{M}_{\odot }/h$ to help decorrelate the parameters a₄, b₄, and c₄.

To ensure that the mass function keeps its approximate power-law behavior when extrapolated to even lower masses (outside our analysis domain), we define a zeroth segment for which we set c₀ = 0 and where a₀ and b₀ are constrained from the first segment in the usual way. The zeroth segment spans the mass range ${10}^{13}\lt M/({M}_{\odot }/h)\lt {10}^{13.1},$ and all subsequent segments span 0.25 dex in mass as discussed.

In practice, for each redshift, we only fit for the c_i that correspond to segments for which the binned halo catalogs do not vanish. All c_i above the highest halo mass are set to the last c_i where binned halo catalogs exist, meaning that the curvature of the log-mass function remains constant toward even larger masses.

3.1.1. Likelihood Function and Halo Mass Function Fits

The input data provided by the simulations for each redshift are the binned halo catalog ${{\boldsymbol{N}}}_{\mathrm{sim}}$ and a covariance matrix ${{\boldsymbol{\Sigma }}}_{\mathrm{sim}}$ (see Section 2.3). As the mass bins are rather coarse, we do not assume the mass function to be constant within a bin, and we explicitly model the number of halos in a given bin i as

$$\begin{eqnarray}&&{N}_{\mathrm{model},i}={L}_{\mathrm{box}}^{3}{\int }_{{M}_{\min ,i}}^{{M}_{\max ,i}}d\mathrm{ln}M\displaystyle \frac{{dn}}{d\mathrm{ln}M},\end{eqnarray} \tag{ 15 }$$

with ${dn}/d\mathrm{ln}M$ from Equation (11). We define the log-likelihood function

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{D}} & \equiv & {{\boldsymbol{N}}}_{\mathrm{sim}}-{{\boldsymbol{N}}}_{\mathrm{model}}\\ \mathrm{ln}{{ \mathcal L }}_{\mathrm{sim}} & = & -\displaystyle \frac{1}{2}{{\boldsymbol{D}}}^{T}{{\boldsymbol{\Sigma }}}_{\mathrm{sim}}^{-1}{\boldsymbol{D}}+\mathrm{const}.\end{array}\end{eqnarray} \tag{ 16 }$$

We regularize the variance in the fit parameters c_i by applying a Gaussian likelihood on the variance between each pair of consecutive c_i,

$$\begin{eqnarray}&&\mathrm{ln}{{ \mathcal L }}_{\mathrm{var}}=-({i}_{\max }-1)\mathrm{ln}\lambda +\displaystyle \sum _{i=1}^{{i}_{\max }-1}-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{c}_{i}-{c}_{i+1}}{\lambda }\right)}^{2},\end{eqnarray} \tag{ 17 }$$

up to a constant. The standard deviation λ is a free parameter of the model. The role of $\mathrm{ln}{{ \mathcal L }}_{\mathrm{var}}$ is to constrain λ and c_i by striking a balance between a good fit (which requires the c_i to differ and thus leads to a noncontinuous second derivative of the log-mass function) and a second derivative of $\mathrm{ln}({dn}/d\mathrm{ln}M)$ that is as little jumpy as possible.

Finally, for each redshift, we constrain the parameters a₄, b₄, c_i, and λ for each design model by sampling the total log-likelihood

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=\mathrm{ln}{{ \mathcal L }}_{\mathrm{sim}}+\mathrm{ln}{{ \mathcal L }}_{\mathrm{var}}+\mathrm{const}.\end{eqnarray} \tag{ 18 }$$

The residuals of the best-fit mass functions obtained this way are consistent with the statistical noise in the halo catalogs. The distribution of reduced χ² across all models and redshifts is close to being centered at unity.

In Figure 3, we show the binned halo catalogs and the mass function fits for model M042 for illustration purposes. The range of best-fit mass functions for all models at z = 0 is shown in the top panel of Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Basis Functions

3.3. (Re)Computing the Principal Component Weights

The PCA performed in the previous subsection provides, for each snapshot, four eigenvectors (or basis functions) and the PC weights associated with each design model. Note that these PC weights are point estimates with no uncertainties, because we performed the PCA on the (zero-uncertainty) best-fit mass function models. To propagate the statistical noise in the halo catalogs onto the PC weights, we now (re)fit the mass function, parameterized by the mean mass function and the four basis functions from the previous subsection, to the halo catalogs. As we are explicitly interested in the uncertainty on the PC weights, we perform a Markov Chain Monte Carlo (MCMC) analysis for each halo catalog using the likelihood function defined in Equation (16).

As a result, the input mass function for every input cosmology and redshift is now described by a four-dimensional posterior parameter distribution. All posterior distributions (8 snapshots × 111 models) are well approximated by multivariate Gaussian distributions, and we thus extract the parameter means and covariances for all models.

We validate the description of the mass functions with these four parameters by studying the residuals between each mean mass function parameterized in this way and its underlying halo catalog. The residuals are shown in Figure 5, which suggests that they are consistent with statistical noise in the halo catalogs. For each redshift, the distribution of residuals across all models is consistent with a reduced χ² of unity.

Zoom In Zoom Out Reset image size

3.4. Gaussian Process Regression

In this section, we briefly review the basics of GP regression and set the appropriate notation. We then apply this method and build an interpolation scheme for the PC weights obtained above.

3.4.1. Basics of Gaussian Process Regression

GP regression assumes a collection of data values ${\boldsymbol{f}}$ at locations ${\boldsymbol{X}}$ that are drawn from a joint Gaussian distribution ${\boldsymbol{f}}={ \mathcal N }(0,{\boldsymbol{K}}({\boldsymbol{X}}))$, with a correlation matrix ${\boldsymbol{K}}$. Without loss of generality, we have subtracted the mean such that the GP is centered at $0$. Following, e.g., Rasmussen & Williams (2006), the GP can be conditioned on training values ${\boldsymbol{f}}$ at inputs ${\boldsymbol{X}}$ with measurement errors ${\boldsymbol{\sigma }}$ such that it predicts function values ${{\boldsymbol{f}}}_{\star }$ at new locations ${{\boldsymbol{X}}}_{\star }$:

$$\begin{eqnarray}&&\begin{array}{l}{\boldsymbol{K}}\equiv {\boldsymbol{K}}({\boldsymbol{X}},{\boldsymbol{X}})\\ {{\boldsymbol{K}}}_{\star }\equiv {\boldsymbol{K}}({{\boldsymbol{X}}}_{\star },{\boldsymbol{X}})\\ {{\boldsymbol{K}}}_{\star \star }\equiv {\boldsymbol{K}}({{\boldsymbol{X}}}_{\star },{{\boldsymbol{X}}}_{\star })\\ {{\boldsymbol{f}}}_{\star }| {{\boldsymbol{X}}}_{\star },{\boldsymbol{X}},{\boldsymbol{f}}\sim \\ { \mathcal N }\left({{\boldsymbol{K}}}_{\star }{\left[{\boldsymbol{K}}+{\boldsymbol{\sigma }}\right]}^{-1}{\boldsymbol{f}},{{\boldsymbol{K}}}_{\star \star }-{{\boldsymbol{K}}}_{\star }{\left[{\boldsymbol{K}}+{\boldsymbol{\sigma }}\right]}^{-1}{{\boldsymbol{K}}}_{\star }\right).\end{array}\end{eqnarray} \tag{ 19 }$$

The correlation matrix ${\boldsymbol{K}}$ is constructed using the correlation function k between two inputs ${\boldsymbol{x}}$ and ${{\boldsymbol{x}}}^{{\prime} }$. A popular choice is the squared exponential kernel, which we adopt in a slightly modified form:

$$\begin{eqnarray}&&k({\boldsymbol{x}},{{\boldsymbol{x}}}^{{\prime} },\lambda ,{\boldsymbol{\rho }})={\lambda }^{-1}\displaystyle \prod _{i=1}^{{N}_{\dim ,{\boldsymbol{x}}}}{\rho }_{i}^{4{\left({x}_{i}-{x}_{i}^{{\prime} }\right)}^{2}},\end{eqnarray} \tag{ 20 }$$

where ${N}_{\dim ,{\boldsymbol{x}}}$ is the dimensionality of ${\boldsymbol{x}}$ and with hyperparameters λ and ${\boldsymbol{\rho }}$, setting the precision of the GP and the spatial correlation length between inputs, respectively. For example, a GP with a large precision λ varies little along spatial directions, resulting in very smooth functions. As ρ_i approaches unity (from below), the correlation length along direction i becomes infinite, and we call the dimension i "inactive"; conversely, for ρ_i ≪ 1, the spatial correlation length becomes very short. The remaining task is thus to set the hyperparameters in an informed way; we will discuss this in Section 3.4.3.

The above summary describes a GP regression scheme for a scalar function f that depends on a multidimensional argument ${\boldsymbol{x}}$. It is straightforward to extend the formalism to vector functions while taking the possible correlations between the function vector elements into account. For ${N}_{\dim ,{\boldsymbol{f}}}$-dimensional data, one constructs ${N}_{\dim ,{\boldsymbol{f}}}$ correlation matrices ${\boldsymbol{K}}$, where each correlation matrix ${{\boldsymbol{K}}}_{n}$ has an independent set of hyperparameters λ_n and ${{\boldsymbol{\rho }}}_{n}$. The measurement errors on the input parameters are now stored in a covariance matrix with ${N}_{\dim ,{\boldsymbol{f}}}$ × ${N}_{\mathrm{design}}$ entries on a side, where ${N}_{\mathrm{design}}$ is the number of input data points. With this formalism, we account for correlated input parameters but assume that the hyperparameters for each input are independent. We will employ such a multidimensional GP regression scheme for our emulator.

3.4.2. Gaussian Process Input Parameters

We prepare the inputs for the GP regression as follows. The cosmological parameters of the design models⁸ are normalized to the range [0, 1]. The mean PC weights from Section 3.3 are standardized to have zero mean and unit variance. The associated covariance matrices are rescaled accordingly.

Thus, we have ${N}_{\dim ,{\boldsymbol{x}}}=8$ function arguments (the cosmological parameters), ${N}_{\dim ,{\boldsymbol{f}}}=4$ dimensional function values (the PC weights), and ${N}_{\mathrm{design}}=111$ input data points.

3.4.3. Hyperparameter Optimization

The parameters ${\lambda }_{1},\ldots ,{\lambda }_{{N}_{\dim ,{\boldsymbol{f}}}}$ and ${{\boldsymbol{\rho }}}_{1},\ldots ,{{\boldsymbol{\rho }}}_{{N}_{\dim ,{\boldsymbol{f}}}}$ are free parameters of the GP model that need to be specified. The marginal likelihood of the GP can be interpreted as a likelihood function of the inputs given the model hyperparameters:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}{ \mathcal L }({\boldsymbol{f}}| {\boldsymbol{X}}) & = & -\displaystyle \frac{1}{2}{{\boldsymbol{f}}}^{T}{\left({\boldsymbol{K}}+{\boldsymbol{\Sigma }}\right)}^{-1}{\boldsymbol{f}}\\ & & -\displaystyle \frac{1}{2}\mathrm{ln}\left[\det ({\boldsymbol{K}}+{\boldsymbol{\Sigma }})\right]+\mathrm{const}.\end{array}\end{eqnarray} \tag{ 21 }$$

Following Habib et al. (2007), we apply priors π on the hyperparameters

$$\begin{eqnarray}\begin{array}{rcl}0.1 & \lt & {\lambda }_{i}\lt 2.5,\\ \pi ({\lambda }_{i}) & \propto & {\lambda }_{i}^{4}\exp (-5{\lambda }_{i}),\quad i=1,\ldots ,{N}_{\dim ,{\boldsymbol{f}}}\end{array}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}\begin{array}{rcl}0.4 & \lt & {\rho }_{{ij}}\lt 1,\\ \pi ({\rho }_{{ij}}) & \propto & {\left(1-{\rho }_{{ij}}\right)}^{-0.9},\quad i=1,\ldots ,{N}_{\dim ,{\boldsymbol{f}}},\\ & & \qquad \qquad \qquad j=1,\ldots ,{N}_{\dim ,{\boldsymbol{x}}}.\end{array}\end{eqnarray} \tag{ 23 }$$

Because we standardized the GP input parameters in the previous subsection, we expect the variance ${\boldsymbol{\lambda }}$ to be close to unity and choose a prior with mean of 1 and standard deviation of 0.45. Similarly, as the design parameters are normalized to [0, 1], the adopted prior on ${\boldsymbol{\rho }}$ encodes our expectation of very smooth functions (and few active dimensions); we have substantial prior mass near 1 and $\Pr ({\rho }_{{ij}}\lt 0.98)\approx 1/3$.

The final step is to set the hyperparameters such that they maximize the likelihood in Equation (21). This represents an inference problem with ${N}_{\dim ,{\boldsymbol{f}}}+{N}_{\dim ,{\boldsymbol{f}}}\times {N}_{\dim ,{\boldsymbol{x}}}=36$ free parameters. For each snapshot, we run an MCMC analysis and set the hyperparameters to the mean recovered values from the MCMC.

3.5. Emulator Output and Uncertainty

For a given requested cosmology, the emulator uses the GP regression scheme to estimate the mean PC weights and their covariance matrix according to Equation (19). As a first consistency check, we confirm that the residual between the four-dimensional PC weights that we obtained by fitting to the halo catalogs of the design models in Section 3.3 is statistically consistent with the PC weights that the GP returns at those locations. This test confirms that the emulator prediction indeed goes through the design data points.

As a next step, we compute the consistency between the PC weights obtained in Section 3.3 for the test models, which were not used for the construction of the emulator, and the prediction by the GP. For three out the eight snapshots, the reduced χ² across the five test models is slightly smaller than unity, indicating that the emulator is making statistically accurate predictions. For the remaining five snapshots, the reduced χ² varies between 1.3 and 2.4. We multiply the GP output covariance matrix with fixed factors such that the reduced χ² is forced to unity. This slightly degrades the emulator precision at the benefit of having statistically consistent residuals with the test models. The emulator is now fully specified.

To obtain predictions for the mass function, the PC weights are multiplied with the mass function basis functions, the mean log-mass function is added, and upon exponentiation, we obtain the emulated mass function (reverse process of Section 3.2). Note that the mass function computed in this way corresponds to the mean predicted PC weights and thus contains no uncertainty. To provide an error estimate, we repeat the computation of the mass function, but instead of taking the mean predicted PC weights, we draw realizations from the multivariate probability distribution that the GP provides, accounting for the additional factors discussed just above. We compute the error on the emulator prediction by taking the variance of such a set of mass function predictions from the GP covariance. Note that the error estimates vary with mass and redshift as they account for the noise in the input halo catalogs; they also vary as a function of the location in cosmology space according to the proximity to design locations (see Figure 6).

Zoom In Zoom Out Reset image size

Note that in principle there is another contribution to the emulator uncertainty due to the uncertainties in the GP hyperparameters. However, across all redshifts, this additional source of noise is about an order of magnitude smaller than the errors discussed above and we therefore neglect it.

3.6. Redshift Evolution of the Halo Mass Function

3.7. Emulator Performance

We verify the emulator performance in two ways: hold-out tests and tests against additional simulations for cosmologies that were not used in the original construction of the emulator.

3.7.1. Hold-out Tests

Hold-out tests are performed by constructing the emulator using all but the ith design model, and this is repeated for each design model i. The hyperparameters of the GP are kept fixed to the values obtained for the full design in Section 3.4.3. Intuitively, the hold-out tests allow us to assess whether the function to be emulated is smooth enough so that ${N}_{\mathrm{design}}-1$ input cosmologies are sufficient for emulation purposes.

The result of the test is shown in Figure 7. The mass functions are correctly predicted at the few-percent level. At high mass, as expected, this test is limited by the statistical noise in the halo catalogs. Note that the hold-out test is powerful, because we are effectively testing the impact of a "hole" in the design, and the full emulator has no such holes. We conclude that our emulator passes this test and that the claim of percent-level mass function predictions made above in Section 3.5 and Figure 6 are justified.

Zoom In Zoom Out Reset image size

3.7.2. Verification with Additional Simulations

An alternative verification scheme consists of carrying out spot checks by comparing the emulator predictions with halo catalogs obtained for additional cosmologies that were not used to construct the emulator; this avoids the "sampling hole" problem (potentially most troublesome near the hypercube boundaries) with hold outs, at the obvious cost of a more systematic coverage of parameter space.

In Figure 8, we compare the emulator prediction with halo catalogs from our five additional cosmologies (M000, T001–T004). Here again, we observe percent-level agreement. By construction, the residuals are within the combined emulator error and the statistical noise in the halo catalogs (see Section 3.5). We thus conclude the verification of the emulator.

Zoom In Zoom Out Reset image size

3.8. Concluding Remarks

We discuss our emulator in the context of existing halo mass function fits and emulators. We consider the simulation for our ΛCDM (and massless-neutrino) cosmology M000, which was not used for the construction of our emulator, and confront it with different predictions from the literature. We perform this comparison at redshift z = 0 and show the results in Figure 9. As discussed above, our emulator provides an accurate prediction for the mass function at the M000 cosmology.

Zoom In Zoom Out Reset image size

A popular set of fitting functions that covers a range of SO mass definitions is described in Tinker et al. (2008); we compare that prediction with our simulation in Figure 9.⁹ Above $M\gtrsim {10}^{14}{M}_{\odot }/h$, this fit provides a very good description of our simulation. At lower mass, however, the Tinker et al. (2008) fit overestimates the number of halos by up to about 10%. This behavior is expected as we extract SO masses for FOF halos, while Tinker et al. (2008) use an SO finder that allows for overlapping subhalos which can boost the number of objects (see the discussion in, e.g., Kravtsov et al. 2004; Tinker et al. 2008; Bocquet et al. 2016). We compare to an alternative fitting function presented in Bocquet et al. (2016) that was calibrated to the Magneticum Pathfinder simulation suite (K. Dolag et al. 2020, in preparation).¹⁰ Importantly, in this work, SO masses were extracted for FOF halos defined by b = 0.16 corresponding to almost the same halo definition as we use (SO masses for b = 0.168 halos). Indeed, the mass function for the Magneticum N-body simulations agrees very well with our M000 model. The mass function for the Magneticum hydrodynamic simulations yields significantly fewer low-mass clusters and groups than their N-body counterparts or our simulation. This is expected, because hydrodynamic feedback effects drive baryons from the halo center, which, for a given halo, leads to a lower total enclosed mass—this is turn leads to a lower halo abundance at fixed mass (e.g., Henson et al. 2017; Springel et al. 2018; see also the discussion in Section 6). In conclusion, after having shown that our emulator passes the verification tests against our own simulations in the previous section, we have now shown that, for our ΛCDM model, our simulation and emulator also agree with two example fitting functions from the literature.¹¹

To date, two emulators for the halo mass function have been presented: the Dark Emulator (Nishimichi et al. 2019) and Aemulus, which is publicly available¹² (McClintock et al. 2019). Both of these are only valid for massless neutrinos and a constant dark energy equation of state while our emulator additionally covers massive neutrinos in the parameter range $0\leqslant {{\rm{\Omega }}}_{\nu }{h}^{2}\leqslant 0.01$ and dynamical dark energy through the w₀–w_a parameterization (see Section 2). As a result, only our emulator can be used for cosmological analyses that aim at constraining the sum of neutrino masses and/or dynamical dark energy.

The Dark Emulator covers a six-dimensional cosmological parameter space. As opposed to our emulator, it is constrained to $\sum {m}_{\nu }=0$ and a constant dark energy equation of state (w_a = 0). Their simulation design contains 101 cosmologies which are sampled from a six-dimensional hypercube. A sliced Latin hypercube design is employed to allow for an even sampling of the parameter space (as a reminder, we achieve even sampling through a tesselation-based nested design strategy; see Section 2). For each model, $1\,\mathrm{Gpc}/h$ and $2\,\mathrm{Gpc}/h$ simulation boxes are run evolving 2048³ particles each; only the smaller boxes are used for the halo mass function emulator. In summary, the Dark Emulator framework is similar to ours except for the more restrictive range of cosmologies.

The Aemulus suite of simulations (DeRose et al. 2019) covers seven cosmological parameters; compared to our work, they impose ${{\rm{\Omega }}}_{\nu }{h}^{2}=0$ and w = const. (equivalent to w_a = 0) but vary the effective number of relativistic species N_eff. Importantly, their design cosmologies are chosen to sample the parameter space allowed by the union of the cosmological constraints from WMAP9+BAO+SN Ia and Planck13+BAO+SN Ia data, and their emulator cannot produce accurate results outside of this region of the parameter space. In contrast, the Mira-Titan Universe design samples the entire parameter hypercube, as described in Section 2. Therefore, this hypercube can be adopted as hard priors in any cosmological analysis that uses our emulator, whereas the Aemulus emulator can only be used in analyses where strong priors from CMB, BAO, and Type Ia supernova measurements are applied. The 40 Aemulus simulations have similar mass resolution as ours but smaller box sizes of $1.05\,\mathrm{Gpc}/h$, thus providing a less precise calibration of the number (density) of higher-mass halos. As a cross-check, we ran the Aemulus emulator for our M000 cosmology. Note that Aemulus provides SO masses ${M}_{200{\rm{m}}}$ (using the rockstar halo finder and thus not the same halo definition as we do); we convert to ${M}_{200{\rm{c}}}$ assuming a Navarro–Frenk–White (NFW) profile (Navarro et al. 1997) and the concentration–mass relation by Child et al. (2018).¹³ As shown in Figure 9 and demonstrated in McClintock et al. (2019), the prediction from Aemulus is very similar to the prediction by Tinker et al. (2008): it agrees with our mass function at about ${10}^{14}{M}_{\odot }/h$ and above (within the noise at high mass) and overestimates the halo abundance toward lower masses by up to about 10%. We attribute the latter effect to different halo definitions, as also discussed in our comparison with Tinker et al. (2008).

The discrepancy between (some of) the mass function predictions below $\sim {10}^{14}{M}_{\odot }/h$, which arises because of different halo (center) definitions, has two important consequences. The first relates to the fact that cosmological analyses of real data must use simulation inputs in a self-consistent manner, i.e., that the mass definition for which the halo mass function is predicted has to be the same mass definition for which the projected halo mass profiles for weak-lensing measurements are used (in other words, all else being equal, the appropriate mass—observable relation will compensate for changes in the mass definition). The second point is more complicated and has to do with halo identification—the halo-finding strategy employed in the simulation must correctly take into account how observed clusters are found in practice (see, e.g., García & Rozo 2019 for additional discussion and an explicit example for optically selected clusters). There must be a consistent pathway to take the halos from a simulated light cone catalog and associate them with single objects that can be found from a simulated observation. (Merging groups and clusters are one example where more care may be needed.) Upcoming and future surveys that push to low cluster masses around ${10}^{14}{M}_{\odot }/h$ and below will need to take this effect into account.

Prior to the emergence of cosmological emulators based on numerical simulations, the cluster cosmology community followed an approach based on the assumed universal form for the mass function. As discussed in the Introduction, in this approach, all of the mass function dependence on cosmology and redshift is fully described as a function of the rms fluctuations σ(M, z) in the linear matter power spectrum $P(k,z,\mathrm{cosmology})$. In this section, we use the Mira-Titan Universe simulations to investigate the limits of the assumed cosmological universality of the mass function for ${M}_{200{\rm{c}}}$.

As a first step, we establish a mass function fit following the universal approach. We choose model M000 (${\rm{\Lambda }}\mathrm{CDM}$ with massless neutrinos) as our baseline model for which we establish the universal fit; we then extrapolate this result to all other Mira-Titan Universe cosmologies and compare to the simulation results.

The universal parameterization of the halo mass function reads

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

with ρ_m designating the matter density. The rms fluctuation σ in the matter density field is related to the linear matter power spectrum P(k, z) via

$$\begin{eqnarray}&&{\sigma }^{2}(M,z)=\displaystyle \frac{1}{2{\pi }^{2}}\int P(k,z){\hat{W}}^{2}({kR}){k}^{2}{dk}\end{eqnarray} \tag{ 25 }$$

with the Fourier transform $\hat{W}$ of the real-space top-hat window function of radius

$$\begin{eqnarray}&&R={\left(\displaystyle \frac{3M}{4\pi {\rho }_{{\rm{m}}}}\right)}^{1/3}.\end{eqnarray} \tag{ 26 }$$

A common parameterization of f(σ) is

$$\begin{eqnarray}&&f(\sigma )=B\left[{\left(\displaystyle \frac{\sigma }{e}\right)}^{-d}+{\sigma }^{-f}\right]\exp (-g/{\sigma }^{2}),\end{eqnarray} \tag{ 27 }$$

with fit parameters d, e, f, and g (e.g., Warren et al. 2006; Tinker et al. 2008) together with a normalization B set by the condition $\int f(\sigma )d\mathrm{ln}{\sigma }^{-1}=1$ which is satisfied for

$$\begin{eqnarray}&&B=2{\left[{e}^{d}{g}^{-d/2}{\rm{\Gamma }}(d/2)+{g}^{-f/2}{\rm{\Gamma }}(f/2)\right]}^{-1}.\end{eqnarray} \tag{ 28 }$$

This model has considerable freedom as it contains an overall amplitude, the slope of the power-law part of the mass function, a relative amplitude of the power-law component, and an exponential cutoff at high mass (low σ).

For cosmological models with massive neutrinos, the degree of universality can be improved by choosing the appropriate power spectrum P(k) that enters Equation (25) and matter density ρ_m that enters Equations (24) and (26). Following, e.g., Ichiki & Takada (2012), Costanzi et al. (2013), and Biswas et al. (2019), we set the matter density to the density of collapsing matter:

$$\begin{eqnarray}&&{\rho }_{{\rm{m}}}={\rho }_{\mathrm{CDM}}+{\rho }_{{\rm{b}}};\end{eqnarray} \tag{ 29 }$$

this is motivated by the fact that neutrinos hardly cluster on halo scales. Following the same reasoning, we use

$$\begin{eqnarray}&&P(k,z)={P}_{\mathrm{CDM}+{\rm{b}}}(k,z){\left(\displaystyle \frac{{\rho }_{\mathrm{CDM}}(z)+{\rho }_{{\rm{b}}}(z)}{{\rho }_{\mathrm{CDM}}(z)+{\rho }_{{\rm{b}}}(z)+{\rho }_{\nu }(z)}\right)}^{2},\end{eqnarray} \tag{ 30 }$$

which is the CDM+baryon power spectrum modified to account for the background density evolution of all particle species.

We obtain parameter posterior distributions for the universal parameters by fitting the mass function from Equation (24) to the M000 halo catalog at z = 0 using the likelihood function from Equation (16). We then use this set of parameters to predict the mass function for every model M001–M111 and the test models T001–T004.

In the top panel of Figure 10, the thick black line indicates that the universal mass function fit obtained for the M000 model indeed provides a good fit to the M000 halo catalog. However, as we use the set of parameters d, e, f, and g to predict the universal mass function for every model M001–M111, large discrepancies become apparent. Up to masses of about ${10}^{14}{M}_{\odot }/h$, the residuals are approximately constant with mass and range from about −20% up to about +30%. At higher masses, as we approach the exponential tail of the mass function, the residuals increase significantly.

Zoom In Zoom Out Reset image size

Note that the degree to which the mass function is universal depends on the specific halo mass definition. For example, the mass function for FOF masses with b = 0.2 or for SO masses ${M}_{180{\rm{m}}}$ was shown to feature a more universal behavior for ${\rm{\Lambda }}\mathrm{CDM}$ cosmologies than the mass function for ${M}_{200{\rm{c}}}$ (White 2002). To extend this statement beyond ${\rm{\Lambda }}\mathrm{CDM}$, Bhattacharya et al. (2011) calibrated an FOF b = 0.2 mass function fit using ${\rm{\Lambda }}\mathrm{CDM}$ simulations and then used a suite of 38 simulations carried out for $w\mathrm{CDM}$ cosmologies to estimate that the universal approach is still valid at the ∼10% level. Additional spot checks using a simulation with a ${\rm{\Lambda }}\mathrm{CDM}$ cosmology with massive neutrinos (νΛCDM) and a simulation with additional dynamical dark energy (νDDE) were performed by Biswas et al. (2019). They found that the universal mass function approach is able to describe the mass function in their νΛCDM cosmology to within 5% and the mass function in their νDDE model at the ∼5% level. In conclusion, we expect that the mass function for ${M}_{200{\rm{c}}}$, which we consider in this work, is less universal than for other mass definitions, and FOF masses with b = 0.2 in particular. Note, however, that the formation and evolution of halos are better understood for SO masses that are defined with respect to the critical density, such as ${M}_{200{\rm{c}}}$ or M_500c, which justifies using these mass definitions despite the mass function being less universal (see, e.g., the discussion in White 2002). Finally, we note that for practical applications, having a sufficiently accurate emulator makes the search for a universal mass function prediction unnecessary. While a universal mass function with a well-characterized uncertainty could be used to expand the cosmological parameter space for which the emulator is valid, the accuracy of the universal mass function outside of the emulator domain of validity would still need to be verified using additional numerical simulations. It remains unclear whether much is gained compared to building a new emulator from such an enhanced simulation campaign.

In the bottom panel of Figure 10, we show the same result as in the upper panel but for the test models T001–T004. As these models were not used for the construction of the emulator, we can directly compare the performance of the universal approach with the emulator. For these models, the universal approach is accurate at the ∼10% level up to about ${10}^{14}{M}_{\odot }/h$ and worsens somewhat toward higher masses. In contrast, the emulator clearly performs much better as its accuracy stays well within 10% up to ${10}^{15}{M}_{\odot }/h$.

We present, and make publicly available,¹⁴ an emulator for the halo mass function for redshifts z ≤ 2 and for masses ${M}_{200{\rm{c}}}\geqslant {10}^{13}{M}_{\odot }/h$ (group and cluster-scale halos) using the newly completed Mira-Titan Universe suite of 111 cosmological N-body simulations. The parameter space spans eight cosmological parameters {${{\rm{\Omega }}}_{{\rm{m}}}{h}^{2}$, ${{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}$, ${{\rm{\Omega }}}_{\nu }{h}^{2}$, ${\sigma }_{8}$, h, n_s, ${w}_{0}$, ${w}_{a}$}. Our work expands the parameter space covered by existing mass function emulators (McClintock et al. 2019; Nishimichi et al. 2019) by also including a dynamical dark energy equation of state and the neutrino mass sum. Following our previous philosophy in terms of sampling parameter space, our simulation design covers the entire parameter hypercube, instead of being centered around the currently available cosmological constraints (as employed, e.g., by the Aemulus emulator). Finally, the use of large simulation boxes (2.1 Gpc) yields better statistics for the high-mass tail of the mass function.

We construct the emulator using GP regression. In contrast to the mass function emulators in the literature, we do not perform the GP regression on the fit parameters of a universal fitting function but on a set of principal component weights. The corresponding PC basis functions are obtained in a data-driven way using piecewise polynomial functions. In our approach, we track the noise in the input halo catalogs and the uncertainties involved in the GP regression. We verify the emulator prediction using hold-out tests and tests against additional simulations that were not used for the construction of the emulator.

The accuracy of the emulator is summarized in Figure 6. The accuracy is a function of mass, redshift, and location in cosmological parameter space. For example, at $M\,={10}^{14}{M}_{\odot }/h$ and z < 1, the error is at the percent level for all cosmologies; at $M={10}^{15}{M}_{\odot }/h$, it stays below 10%.

As pointed out in Section 4, it is important that observational studies use consistent halo definitions for all simulation-based data products that enter the analysis. Furthermore, studies of deep cluster surveys that extend down to ${10}^{14}{M}_{\odot }/h$ or below need to ensure that the halos identified in the simulations match the objects identified in the data. For example, the halo finder may identify two objects (that may partially overlap) when a survey with finite spatial resolution would discover a single, merged cluster. Forward modeling the cluster abundance as a function of the observable quantities (instead of going through the halo mass function) directly from the simulations may be a promising way to address the issue.

This work does not address the influence of nongravitational feedback effects on the mass function. This feedback, which is mainly driven by active galactic nuclei, redistributes the material in the cluster center and thereby alters the halo mass profile (e.g., Henson et al. 2017; Springel et al. 2018). The impact on low-mass clusters and groups is more important than on massive clusters, which results in the mass function being much less affected in the high-mass regime (e.g., Cusworth et al. 2014; Schaller et al. 2015; Bocquet et al. 2016). As a first step, our gravity-only mass function may be modified in postprocessing to approximately account for the impact of baryonic effects. Recent work using a set of six hydrodynamic BAHAMAS simulations, which include the effects of dynamical dark energy, suggests that the impact of baryonic and cosmological effects on the mass function can be separated with an accuracy of ≈1%–2% (Pfeifer et al. 2020). It remains to be seen whether such a modified mass function prediction is accurate enough for cosmological analyses of future, high-quality data sets. If not, the mass function prediction will have to be built from large suites of hydrodynamic simulations, which are expected to be available in the near future. In this case, the systematic uncertainty in the hydrodynamic modeling will have to be accounted for and described along with the cosmological parameters in a rigorous statistical analysis framework. The emulator framework presented in this work represents an important step in this direction.

We thank Derek Bingham for many discussions and for his contributions to the tessellation-based sampling scheme used here. We also record our debt to Dave Higdon for generously sharing many contributions and sharp insights for more than a decade. We further thank Jörg Dietrich and Sebastian Grandis for fruitful discussions about vector-function Gaussian processes and the mass function fits, respectively. We thank Takahiro Nishimichi and the Dark Emulator team for sharing some of their input mass functions for cross-checking purposes. S.B. acknowledges the support by the ORIGINS Cluster (funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy—EXC-2094—390783311) and the Ludwig-Maximilians-Universität München. S.H., K.H., and E.L. acknowledge support from the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, and Office of High Energy Physics, Scientific Discovery through Advanced Computing (SciDAC) program under award No. 231018. The work of the authors at Argonne National Laboratory was supported under the U.S. Department of Energy contract DE-AC02-06CH11357. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC02-06CH11357. This research also used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC05-00OR22725.

Software: astropy (Astropy Collaboration et al. 2018), GetDist,¹⁵ numpy (van der Walt et al. 2011), scipy, matplotlib (Hunter 2007), polychord (Handley et al. 2015), pyGTC (Bocquet & Carter 2016), (py)MultiNest (Buchner et al. 2014; Feroz et al. 2019).

We present all cosmologies of the Mira-Titan Universe and of the additional test models in Table A1.

Table A1.  Cosmologies in the Mira-Titan Universe

Model	${{\rm{\Omega }}}_{{\rm{m}}}{h}^{2}$	${{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}$	σ8	h	ns	w0	wa	${{\rm{\Omega }}}_{\nu }{h}^{2}$
Massless-neutrino design
M001	0.1472	0.02261	0.8778	0.6167	0.9611	−0.7000	0.67220	0.0
M002	0.1356	0.02328	0.8556	0.7500	1.0500	−1.0330	0.91110	0.0
M003	0.1550	0.02194	0.9000	0.7167	0.8944	−1.1000	−0.28330	0.0
M004	0.1239	0.02283	0.7889	0.5833	0.8722	−1.1670	1.15000	0.0
M005	0.1433	0.02350	0.7667	0.8500	0.9833	−1.2330	−0.04445	0.0
M006	0.1317	0.02150	0.8333	0.5500	0.9167	−0.7667	0.19440	0.0
M007	0.1511	0.02217	0.8111	0.8167	1.0280	−0.8333	−1.00000	0.0
M008	0.1200	0.02306	0.7000	0.6833	1.0060	−0.9000	0.43330	0.0
M009	0.1394	0.02172	0.7444	0.6500	0.8500	−0.9667	−0.76110	0.0
M010	0.1278	0.02239	0.7222	0.7833	0.9389	−1.3000	−0.52220	0.0
Main design
M011	0.1227	0.0220	0.7151	0.5827	0.9357	−1.0821	1.0646	0.000345
M012	0.1241	0.0224	0.7472	0.8315	0.8865	−1.2325	−0.7646	0.001204
M013	0.1534	0.0232	0.8098	0.7398	0.8706	−1.2993	1.2236	0.003770
M014	0.1215	0.0215	0.8742	0.5894	1.0151	−0.7281	−0.2088	0.001752
M015	0.1250	0.0224	0.8881	0.6840	0.8638	−1.0134	0.0415	0.002789
M016	0.1499	0.0223	0.7959	0.6452	1.0219	−1.0139	0.9434	0.002734
M017	0.1206	0.0215	0.7332	0.7370	1.0377	−0.9472	−0.9897	0.000168
M018	0.1544	0.0217	0.7982	0.6489	0.9026	−0.7091	0.6409	0.006419
M019	0.1256	0.0222	0.8547	0.8251	1.0265	−0.9813	−0.3393	0.004673
M020	0.1514	0.0225	0.7561	0.6827	0.9913	−1.0101	−0.7778	0.009777
M021	0.1472	0.0221	0.8475	0.6583	0.9613	−0.9111	−1.5470	0.000672
M022	0.1384	0.0231	0.8328	0.8234	0.9739	−0.9312	0.5939	0.008239
M023	0.1334	0.0225	0.7113	0.7352	0.9851	−0.8971	0.3247	0.003733
M024	0.1508	0.0229	0.7002	0.7935	0.8685	−1.0322	1.0220	0.003063
M025	0.1203	0.0230	0.8773	0.6240	0.9279	−0.8282	−1.5005	0.007024
M026	0.1224	0.0222	0.7785	0.7377	0.8618	−0.7463	0.3647	0.002082
M027	0.1229	0.0234	0.8976	0.8222	0.9698	−1.0853	0.8683	0.002902
M028	0.1229	0.0231	0.8257	0.6109	0.9885	−0.9311	0.8693	0.009086
M029	0.1274	0.0228	0.8999	0.8259	0.8505	−0.7805	0.5688	0.006588
M030	0.1404	0.0222	0.8232	0.6852	0.8679	−0.8594	−0.4637	0.008126
M031	0.1386	0.0229	0.7693	0.6684	1.0478	−1.2670	1.2536	0.006502
M032	0.1369	0.0215	0.8812	0.8019	1.0005	−0.7282	−1.6927	0.000905
M033	0.1286	0.0230	0.7005	0.6752	1.0492	−0.7119	−0.8184	0.007968
M034	0.1354	0.0216	0.7018	0.5970	0.8791	−0.8252	−1.1148	0.003620
M035	0.1359	0.0228	0.8210	0.6815	0.9872	−1.1642	−0.1801	0.004440
M036	0.1390	0.0220	0.8631	0.6477	0.8985	−0.8632	0.8285	0.001082
M037	0.1539	0.0224	0.8529	0.5965	0.8943	−1.2542	0.8868	0.003549
M038	0.1467	0.0227	0.7325	0.5902	0.9562	−0.8019	0.3628	0.007077
M039	0.1209	0.0223	0.8311	0.7327	0.9914	−0.7731	0.4896	0.001973
M040	0.1466	0.0229	0.8044	0.8015	0.9376	−0.9561	−0.0359	0.000893
M041	0.1274	0.0218	0.7386	0.6752	0.9707	−1.2903	1.0416	0.003045
M042	0.1244	0.0230	0.7731	0.6159	0.8588	−0.9043	0.8095	0.009194
M043	0.1508	0.0233	0.7130	0.8259	0.9676	−1.0551	0.3926	0.009998
M044	0.1389	0.0224	0.8758	0.6801	0.9976	−0.8861	−0.1804	0.008018
M045	0.1401	0.0228	0.7167	0.6734	0.9182	−1.2402	1.2155	0.006610
M046	0.1381	0.0224	0.7349	0.8277	1.0202	−1.1052	−1.0533	0.006433
M047	0.1411	0.0216	0.7770	0.7939	0.9315	−0.8042	0.7010	0.003075
M048	0.1374	0.0226	0.7683	0.6865	0.8576	−1.1374	−0.5106	0.004548
M049	0.1339	0.0217	0.7544	0.5920	1.0088	−0.8520	−0.7438	0.003512
M050	0.1337	0.0233	0.8092	0.7309	0.9389	−0.7230	0.6920	0.005539
M051	0.1514	0.0222	0.7433	0.6502	0.8922	−0.9871	0.8803	0.002842
M052	0.1483	0.0230	0.7012	0.6840	0.9809	−1.2881	−0.9045	0.006199
M053	0.1226	0.0226	0.7998	0.8265	1.0161	−1.2593	−0.3858	0.001096
M054	0.1345	0.0216	0.8505	0.6251	0.8535	−1.2526	0.5703	0.007438
M055	0.1298	0.0222	0.7504	0.8170	0.9574	−1.0573	1.0338	0.006843
M056	0.1529	0.0219	0.8508	0.6438	1.0322	−0.7359	0.6931	0.006311
M057	0.1419	0.0234	0.7937	0.7415	1.0016	−0.7710	−1.5964	0.005128
M058	0.1226	0.0224	0.7278	0.6152	1.0348	−1.1051	0.2955	0.007280
M059	0.1529	0.0224	0.7035	0.6877	0.8616	−0.9833	−1.1788	0.009885
M060	0.1270	0.0233	0.8827	0.5622	0.8609	−1.1714	0.8346	0.009901
M061	0.1272	0.0220	0.8021	0.8302	0.8968	−0.9545	−0.6659	0.004782
M062	0.1257	0.0216	0.7699	0.5813	0.9460	−0.8041	0.7956	0.003922
M063	0.1312	0.0229	0.7974	0.5890	0.9522	−0.9560	0.6650	0.001740
M064	0.1437	0.0218	0.8866	0.7402	0.9335	−1.0713	0.7128	0.003782
M065	0.1445	0.0218	0.8797	0.5554	0.9353	−1.2423	−1.3032	0.008850
M066	0.1392	0.0229	0.8849	0.8176	0.8579	−1.2334	0.7098	0.000187
M067	0.1398	0.0224	0.7795	0.7473	1.0392	−1.0471	0.7377	0.008256
M068	0.1319	0.0232	0.7645	0.8112	0.9199	−0.7812	0.1646	0.003716
M069	0.1392	0.0227	0.8130	0.6062	0.8765	−1.0792	0.8817	0.006372
M070	0.1499	0.0232	0.8093	0.7587	0.9260	−0.8942	−0.9967	0.009760
M071	0.1218	0.0224	0.7348	0.7999	1.0330	−0.9341	0.8896	0.002212
M072	0.1265	0.0218	0.8349	0.6080	0.9291	−1.1293	0.2197	0.002806
M073	0.1212	0.0227	0.7683	0.6587	0.8704	−0.9662	0.9453	0.000327
M074	0.1337	0.0216	0.8575	0.8099	0.8518	−1.0815	1.0506	0.002370
M075	0.1484	0.0230	0.8491	0.7212	0.9566	−0.8980	0.8181	0.002716
M076	0.1419	0.0215	0.7700	0.6091	0.9332	−0.9753	−0.2691	0.008143
M077	0.1321	0.0230	0.7011	0.6562	0.9938	−1.1170	1.0706	0.001979
M078	0.1446	0.0219	0.7903	0.8074	0.9752	−1.2323	1.1681	0.004021
M079	0.1344	0.0235	0.8742	0.7575	0.9219	−1.0483	−0.3793	0.004423
M080	0.1419	0.0218	0.8420	0.8205	0.9145	−1.1295	−1.2357	0.001959
M081	0.1366	0.0224	0.7034	0.6599	0.8745	−0.8122	0.7675	0.005665
M082	0.1511	0.0218	0.8061	0.7242	1.0132	−0.7940	0.0740	0.004488
M083	0.1501	0.0230	0.7458	0.6037	0.9999	−1.2300	0.9126	0.008023
M084	0.1244	0.0227	0.8444	0.7462	1.0351	−1.2012	1.0042	0.002919
M085	0.1546	0.0227	0.8201	0.8187	0.8620	−1.0794	0.3306	0.005524
M086	0.1289	0.0221	0.8467	0.7498	0.9158	−0.8964	0.7044	0.006605
M087	0.1437	0.0235	0.8383	0.6612	1.0206	−0.7128	0.3537	0.006952
M088	0.1237	0.0229	0.8779	0.6050	0.8725	−1.2333	1.1145	0.001034
M089	0.1505	0.0221	0.8396	0.5830	0.8506	−0.8262	0.3234	0.002603
M090	0.1439	0.0226	0.8366	0.6987	0.9116	−1.2874	0.2509	0.009071
M091	0.1325	0.0224	0.8076	0.6680	0.9435	−0.7361	−0.9543	0.003494
M092	0.1476	0.0215	0.8452	0.8060	0.9855	−0.9543	0.9158	0.007598
M093	0.1389	0.0219	0.7151	0.6105	0.9228	−1.2533	−0.3018	0.004566
M094	0.1386	0.0235	0.7699	0.7494	0.8529	−1.1243	1.0953	0.006593
M095	0.1424	0.0223	0.8764	0.6612	0.9422	−1.2912	1.2729	0.002028
M096	0.1404	0.0219	0.8946	0.8155	1.0442	−1.1562	−0.8081	0.001851
M097	0.1351	0.0226	0.7560	0.6549	1.0041	−0.8389	0.8124	0.005556
M098	0.1259	0.0225	0.7918	0.7512	0.9055	−1.1744	0.9019	0.003027
M099	0.1237	0.0233	0.8907	0.6375	0.9715	−1.2563	−0.6293	0.007970
M100	0.1232	0.0221	0.7715	0.5530	0.8635	−0.9174	−1.7275	0.007150
M101	0.1526	0.0217	0.7535	0.7292	0.8836	−0.7672	−0.1455	0.004596
M102	0.1531	0.0229	0.8727	0.8137	0.9916	−1.1062	0.5227	0.005416
M103	0.1274	0.0223	0.8993	0.7448	1.0455	−0.9231	0.7886	0.006497
M104	0.1481	0.0222	0.7512	0.7255	1.0029	−1.0720	0.1433	0.000910
M105	0.1490	0.0222	0.8922	0.5780	0.9802	−0.8530	0.4716	0.002495
M106	0.1422	0.0223	0.8679	0.6048	0.8869	−0.8012	0.6662	0.009949
M107	0.1304	0.0233	0.8171	0.8062	1.0495	−0.8080	0.3337	0.003607
M108	0.1414	0.0223	0.7269	0.7524	0.9095	−1.0203	0.6065	0.008364
M109	0.1377	0.0229	0.8656	0.6012	1.0062	−1.1060	0.9672	0.006263
M110	0.1336	0.0220	0.8703	0.8423	0.9509	−1.1045	0.0875	0.009305
M111	0.1212	0.0230	0.7810	0.6912	0.9695	−0.9892	0.1224	0.007263
Additional test models
T001	0.1333	0.02170	0.8233	0.7444	0.9778	−1.1560	−1.1220	0.005311
T002	0.1450	0.02184	0.8078	0.6689	0.9000	−0.9333	−0.5667	0.003467
T003	0.1417	0.02300	0.7767	0.7256	0.9222	−0.8444	0.8222	0.004389
T004	0.1317	0.02242	0.7611	0.6878	0.9333	−1.2000	−0.2889	0.001622
M000	0.1335	0.02258	0.8	0.71	0.963	−1.0	0.0	0.0

Note. The models are chosen according to a tesselation-based nested design strategy (Heitmann et al. 2016). The models M001–M010 form a complete design with massless neutrinos. Within the main design, the first 26 models (M011–M036), the first 55 models (M011–M066), and the full design (M011–M111), form complete designs, with improving refinement. The emulator presented in this work is constructed using M001–M111.

Download table as:  ASCIITypeset images: 1 2