The Mira-Titan Universe. III. Emulation of the Halo Mass Function
Sebastian Bocquet Katrin Heitmann Salman Habib Earl Lawrence Thomas Uram Nicholas Frontiere Adrian Pope Hal Finkel
astro-ph.CO





0000-0002-4900-805X]Sebastian Bocquet
sebastian.bocquet@physik.lmu.de



0000-0003-1468-8232]Katrin Heitmann

0000-0002-7832-0771]Salman Habib













ALCF Division, Argonne National Laboratory, Lemont, IL 60439, USA
CPS Division, Argonne National Laboratory, Lemont, IL 60439, USA
HEP Division, Argonne National Laboratory, Lemont, IL 60439, USA
CCS-6, CCS Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 Munich, Germany
Excellence Cluster ORIGINS, Boltzmannstr. 2, 85748 Garching, Germany




We construct an emulator for the halo mass function over group and cluster mass scales for a range of cosmologies, including the effects of dynamical dark energy and massive neutrinos. The emulator is based on the recently completed Mira-Titan Universe suite of cosmological N-body simulations. The main set of simulations spans 111 cosmological models with 2.1 Gpc boxes. We extract halo catalogs in the redshift range z=[0.0, 2.0] and for masses M_200c≥ 10^13M_⊙/h. The emulator covers an 8-dimensional hypercube spanned by {Ω_mh^2, Ω_bh^2, Ω_ν h^2, σ_8, h, n_s, w_0, w_a}; spatial flatness is assumed. We obtain smooth halo mass functions by fitting piecewise second-order polynomials to the halo catalogs and employ Gaussian process regression to construct the emulator while keeping track of the statistical noise in the input halo catalogs and uncertainties in the regression process. For redshifts z≲1, the typical emulator precision is better than 2% for 10^13-10^14 M_⊙/h and <10% for M≃ 10^15M_⊙/h. For comparison, fitting functions using the traditional universal form for the halo mass function can be biased at up to 30% at M≃ 10^14M_⊙/h for z=0. Our emulator is publicly available at <https://github.com/SebastianBocquet/MiraTitanHMFemulator>.







§ INTRODUCTION


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 <cit.>. 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 (; for a review, see ). In practice, such halos are detected via the (localized) presence of galaxies (, ), gravitational lensing signatures <cit.>, 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 (, , ) and Sunyaev-Zel'dovich effect observations (, , )).
The cosmological constraints resulting from the analyses of cluster samples already allow for competitive measurements of the growth history of the universe <cit.>. 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 <cit.>.

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., ). 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 <cit.>, further formalized via the excursion set approach <cit.>. 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 <cit.>. 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 <cit.>. (For a discussion on the possible halo definition-dependence of the universal form, see .)

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 <cit.>.
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 <cit.>.

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 <cit.>, especially if one wishes to include the effects due to non-zero 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 <cit.>. 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 <cit.>.

In this paper, we present a mass function emulator using the  suite of N-body simulations <cit.>. 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^3 particles) of the simulations were designed to measure the mass function at mass scales starting at ∼10^13 M_⊙/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^15M_⊙/h. This new emulator is a complement to the Mira-Titan emulator for the matter power spectrum presented in <cit.>.

This paper is structured as follows. The  simulations and the extraction of binned halo catalogs are discussed in Section <ref>. We describe the construction of the halo mass function emulator and verify its accuracy in Section <ref>. We discuss existing predictions for the mass function in Section <ref>. In Section <ref>, we use the  simulations to discuss the limits of cosmological universality of the mass function. We conclude with a summary and outlook in Section <ref>.






§ THE MIRA-TITAN UNIVERSE



 The simulations used in this work were run using the HACC cosmological N-body code <cit.>. The suite of simulations is named the  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 <cit.>. Here, we summarize the relevant aspects, and refer the reader to the original paper and references therein for further details.



 §.§ Cosmological Modeling: ν w_0 w_aCDM


The  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 {, , , , h, n_s, , }. All models are spatially flat (Ω_k=0).
We consider a dynamical dark energy equation of state, using the common parameterization <cit.>

    w(a) = w_0 + w_a(1-a).

We further define the parameter combination

    w_b ≡ (-w_0-w_a)^1/4

as we will not consider the (w_0,w_a) parameter space but (w_0,w_b) (see discussion after Eq. <ref>).

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 <cit.>. A detailed description of our approach is given in <cit.> and <cit.>. Its validity with regard to power spectrum measurements on large to quasi-linear scales is discussed extensively in <cit.> and the effect of neutrinos on the mass function is investigated for two models in <cit.>, 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, the neutrino clustering itself is not taken into account. Particular care is given to the set-up 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 σ_8. 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 since 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 σ_8. 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 neglected. 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 sub-dominant 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 post-processing gravity-only simulations) will be added over time <cit.>.



lccccc

 simulations used to construct the emulator.
Model     Box size     N_particle     Force res.
    M_particle     min()
M001–M111     2.100 Gpc     3200^3     6.6 kpc     1.08×10^10(Ω_m - Ω_ν/0.28)(h/0.7)^2M_⊙     10^13M_⊙/h 
M006     2.091 Gpc     3200^3     6.6 kpc     1.02×10^10M_⊙     10^13M_⊙/h 

M023     2.085 Gpc     3200^3     6.6 kpc     9.96×10^9M_⊙     10^13M_⊙/h 

M046     1.865 Gpc     3200^3     6.6 kpc     7.23×10^9M_⊙     10^13M_⊙/h
The mass cuts min() 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 this change in volume is trivial to account for.




 §.§ The Design


The  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  has an inbuilt notion of sequential convergence and is aimed at generating multiple emulators at percent-level error levels, including the halo mass function <cit.>. 





The cosmological parameters are chosen within the following boundaries:

    0.12≤   Ω_m h^2   ≤ 0.155
    
    0.0215≤   Ω_b h^2   ≤ 0.0235
    
    0.0≤   Ω_ν h^2   ≤ 0.01
    
    0.7≤   σ_8   ≤ 0.9
    
    0.55≤   h   ≤ 0.85
    
    0.85≤   n_s   ≤ 1.05
    
    -1.3≤   w_0   ≤ -0.7
    
    0.3 ≤   w_b   ≤ 1.3

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 <cit.>. As a result, w_a is jointly constrained with w_0; the smallest allowed value is w_a=-2.16 and the highest allowed value is w_a=1.29.
The range in  corresponds to a range in the sum of neutrino masses 0≤≤0.94 eV. While the upper limit is significantly higher than current results for a  background cosmology, the constraint on  significantly degrades when w_0 and w_a are also allowed to vary (see Figure <ref>).

The choice of the 111 design models and their space-filling properties are discussed in Section 3 in <cit.>. 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 =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 <ref> and are visualized in Figure <ref>.[Figure <ref> also shows current cosmological constraints from Planck TTTEEE+lowl+lowE <cit.>, Baryon Acoustic Oscillations <cit.>, and the Pantheon supernova sample <cit.>.] As discussed earlier, the design hypercube is constructed using w_0 and w_b, which leads to the impression of uneven sampling in (w_0,w_a) space.

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

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 were 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 <cit.>. All of the five reference simulations evolve 3200^3 particles as in our main simulation suite.







 §.§ 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 8 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^13M_⊙/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 Δlog_10M = 0.1. Then, we combine the highest-mass bins such that the resulting (wide) high-mass bin contains at least 20 halos.
Using 4^3 spatial jackknife sub-volumes, 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 <ref>, we show the binned halo mass functions for the 111 design cosmologies.






§ EMULATOR CONSTRUCTION



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

    
  * For every input cosmology, obtain a smooth halo mass function model from the halo catalog 
    
  * Perform a principal component analysis (PCA) on the (logarithm of the) smooth mass function models
    
  * Keep the first four eigenvectors as basis functions
    
  * Fit the four basis functions to each halo catalog to obtain a four-dimensional posterior parameter distribution for each model
    
  * Set up Gaussian Process regression to interpolate between the four fit parameters, accounting for their covariances
    
  * Set the hyperparameters of the Gaussian Process
    
  * 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_⊙/h and volume in units of (Mpc/h)^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.






 §.§ 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

    ln(dn/dln M)_i = a_i + b_i ln(M/M_piv) + c_i ln(M/M_piv)^2.

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

    a_i    + b_i ln(M_join/M_piv) + c_i ln(M_join/M_piv)^2 = 
       a_i+1 + b_i+1ln(M_join/M_piv) + c_i+1ln(M_join/M_piv)^2

and its derivative

    b_i + 2c_i ln(M_join/M_piv) = b_i+1 + 2c_i+1ln(M_join/M_piv).

In addition, we require 

    c_i<0

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_4 as well as b_4 and c_4. 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 Eqs. <ref> and <ref>. We define M_piv=10^14M_⊙/h to help de-correlate the parameters a_4, b_4, and c_4.

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 0-th segment for which we set c_0=0 and where a_0 and b_0 are constrained from the first segment in the usual way. The 0-th segment spans the mass range 10^13<M/(M_⊙/h)<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.





  §.§.§ Likelihood Function and Halo Mass Function Fits


The input data provided by the simulations for each redshift are the binned halo catalog N⃗_sim and a covariance matrix Σ⃗_sim (see Section <ref>). 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

    N_model, i = L_box^3 ∫_M_min, i^M_max, idln M dn/dln M

with dn/dln M from Eq. <ref>. We define the log-likelihood function

    D⃗   ≡N⃗_sim - N⃗_model
    lnℒ_sim   = -1/2D⃗^T Σ⃗_sim^-1D⃗ + const.


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

    lnℒ_var = -(i_max-1)lnλ + ∑_i=1^i_max-1 -1/2(c_i-c_i+1/λ)^2

up to a constant. The standard deviation λ is a free parameter of the model. The role of lnℒ_var is to constrain λ and the c_i by striking a balance between a good fit (which requires the c_i to differ and thus leads to a non-continuous second derivative of the log-mass function) and a second derivative of ln(dn/dln M) that is as little jumpy as possible.

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

    lnℒ = lnℒ_sim + lnℒ_var + const.


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 χ^2 across all models and redshifts is close to being centered at unity.

In Figure <ref>, 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 <ref>.










 §.§ Basis Functions



From the previous subsection, we now have 111 best-fit mass function models for each snapshot. To compress that information, we perform a PCA. We prepare the best-fit mass function models by taking their logarithm to reduce the dynamic range and by subtracting the mean log-model (see second panel of Figure <ref>). Because of the exponential high-mass cutoff of the mass function, it is clear that the best-fit mass functions diverge in the high-mass regime. As we do not want this divergence to dominate the PC decomposition, we limit the best-fit mass functions to the mass range covered by the binned halo catalogs. After carrying out some sensitivity tests, we decided to keep the first four principal components and their associated four eigenvectors in mass function space. With four components, the error introduced by the PCA is about 1% – although this value is not directly relevant as we re-compute all PC weights in the next subsection. The eigenvectors and eigenvalues for z=0 are shown in the two bottom panels of Figure <ref>.






 §.§ (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, parametrized 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 Eq. <ref>.

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 multi-variate 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 parametrized in this way and its underlying halo catalog. The residuals are shown in Figure <ref> 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 χ^2 of unity.







 §.§ Gaussian Process Regression

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






  §.§.§ Basics of Gaussian Process Regression


Gaussian Process regression assumes a collection of data values f⃗ at locations X⃗ that are drawn from a joint Gaussian distribution f⃗ = 𝒩(0⃗,K⃗(X⃗)), with a correlation matrix K⃗. Without loss of generality, we have subtracted the mean such that the GP is centered at 0⃗. Following, e.g., <cit.>, the GP can be conditioned on training values f⃗ at inputs X⃗ with measurement errors σ⃗ such that it predicts function values f⃗_⋆ at new locations X⃗_⋆:

    K⃗≡K⃗(X⃗,X⃗)
       K⃗_⋆≡K⃗(X⃗_⋆,X⃗)
       K⃗_⋆⋆≡K⃗(X⃗_⋆,X⃗_⋆)
       f⃗_⋆|X⃗_⋆, X⃗, f⃗∼
       𝒩(K⃗_⋆[K⃗ + σ⃗]^-1f⃗,  K⃗_⋆⋆ - K⃗_⋆[K⃗ + σ⃗]^-1K⃗_⋆).


The correlation matrix K⃗ is constructed using the correlation function k between two inputs x⃗ and x⃗^'. A popular choice is the squared exponential kernel, which we adopt in a slightly modified form

    k(x⃗, x⃗^', λ, ρ⃗) = λ^-1∏_i=1^ρ_i^4(x_i-x_i^')^2,

where  is the dimensionality of x⃗ and
with hyperparameters λ and ρ⃗, 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 <ref>.

The above summary describes a GP regression scheme for a scalar function f that depends on a multi-dimensional argument x⃗. It is straightforward to extend the formalism to vector functions while taking the possible correlations between the function vector elements into account. For -dimensional data, one constructs  correlation matrices K⃗, where each correlation matrix K⃗_n has an independent set of hyperparameters λ_n and ρ⃗_n. The measurement errors on the input parameters are now stored in a covariance matrix with × entries on a side, where  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 multi-dimensional GP regression scheme for our emulator.






  §.§.§ GP Input Parameters

We prepare the inputs for the GP regression as follows. The cosmological parameters of the design models[We choose the w_0-w_b parametrization of the dark energy equation of state (see Eq. <ref>).] are normalized to the range [0,1]. The mean PC weights from Section <ref> are standardized to have zero mean and unit variance. The associated covariance matrices are rescaled accordingly.

Thus, we have =8 function arguments (the cosmological parameters), =4 dimensional function values (the PC weights), and =111 input data points.






  §.§.§ Hyperparameter Optimization


The parameters λ_1,…,λ_ and ρ⃗_1,…,ρ⃗_ 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:

    lnℒ(f⃗|X⃗) =    -1/2f⃗^T (K⃗ + Σ⃗)^-1f⃗
       - 1/2ln[(K⃗ + Σ⃗)] + const.

Following <cit.>, we apply priors π on the hyperparameters

    0.1    < λ_i < 2.5, 
    π(λ_i)    ∝λ_i^4 exp(-5λ_i),     i=1,…,
    
        0.4    < ρ_ij < 1, 
    π(ρ_ij)    ∝ (1-ρ_ij)^-0.9,     i=1,…,, 
          j=1,…,.

Because we standardized the GP input parameters in the previous subsection, we expect the variance λ⃗ 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 ρ⃗ encodes our expectation of very smooth functions (and few active dimensions); we have substantial prior mass near 1, and Pr(ρ_ij<0.98)≈1/3.

The final step is to set the hyperparameters such that they maximize the likelihood in Eq. <ref>. This represents an inference problem with +×=36 free parameters. For each snapshot, we run an MCMC analysis and set the hyperparameters to the mean recovered values from the MCMC.







 §.§ 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 Eq. <ref>. 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 <ref> are 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 <ref> 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 χ^2 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 χ^2 vary between 1.3 and 2.4. We multiply the GP output covariance matrix with fixed factors such that the reduced χ^2 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 <ref>). 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 <ref>).



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.







 §.§ Redshift Evolution of the Halo Mass Function

Our emulator provides the mass function for the requested cosmology for 8 discrete redshifts between 0 and 2. To obtain the mass function at intermediate redshifts, we recommend linearly interpolating the logarithm of the mass function ln[dn(M,z)/dln M].





 §.§ 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.









  §.§.§ Hold-Out Tests


Hold-out tests are performed by constructing the emulator using all but the i-th 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 <ref>. Intuitively, the hold-out tests allow us to assess whether the function to be emulated is smooth enough so that N_design-1 input cosmologies are sufficient for emulation purposes.

The result of the test is shown in Figure <ref>. 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 <ref> and Figure <ref> are justified.








  §.§.§ 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 <ref>, 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 <ref>). We thus conclude the verification of the emulator.








 §.§ Concluding Remarks


We conclude this section with a discussion of the approach adopted for the construction of the emulator. The initial mass function fits described in Section <ref> are only used to construct a (close to) optimal set of basis functions in the subsequent Section <ref>. We then keep four basis functions and obtain four-dimensional parameter sets by fitting these basis functions to the halo catalogs in Section <ref>. In principle, any other set of basis functions would work, too. One could, for example, directly use the binned halo catalogs as inputs to the PCA. This would, however, lead to non-continuous mass function predictions and higher overall noise levels which is why we did not follow this approach.

In yet another possible analysis setup, we used the universal mass function functional form (described in detail in Section <ref>) to obtain smooth mass function fits for each model. However, we found that this parametrized mass function is not able to accurately capture the details of the behavior at high mass. The inaccuracy is very subtle: for a given model, the residuals and reduced χ^2 are completely acceptable. However, when measuring the residuals for all 111 models, we noticed that the last one or two bins in the halo catalogs almost all scattered high compared to their respective best-fit mass function model. The distribution of reduced χ^2 for all models indeed suggested a slightly bad fit.

By adopting the approach of using a constrained piecewise polynomial model instead of raw binned halo counts or the universal fitting function approach, we strike a balance between being purely data-driven and incorporating physical intuition about the mass function. In summary, our mass function emulator is built upon the assumptions that the (logarithm of the) mass function can be described by a second-order polynomial in log-mass whose second derivative varies on scales of about 0.25 dex in mass, and that the variation of the (re-fitted) PC weights with cosmology can be captured by a multi-variate Gaussian distribution. The residuals with respect to the design halo catalogs were shown in Figure <ref>.






§ DISCUSSION OF EXISTING HALO MASS FUNCTION FITS AND EMULATORS



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 <ref>. As discussed above, our emulator provides an accurate prediction for the mass function at the M000 cosmology.

A popular set of fitting functions that covers a range of SO mass definitions is described in <cit.>; we compare that prediction with our simulation in Figure <ref>.[We verify our implementation of the fitting function against the HMFcalc tool and find agreement better than 0.4% <cit.>.] Above M≳10^14M_⊙/h, this fit provides a very good description of our simulation. At lower mass, however, the <cit.> fit overestimates the number of halos by up to about 10%. This behavior is expected since we extract SO masses for FOF halos, while <cit.> use an SO finder that allows for overlapping subhalos which can boost the number of objects <cit.>. We compare to an alternative fitting function presented in <cit.> that was calibrated to the Magneticum Pathfinder simulation suite (Dolag et al., in prep.).[<http://www.magneticum.org/>] 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 less low-mass clusters and groups than their N-body counterparts or our simulation. This is expected, since 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 <cit.>. 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.[In Section <ref>, we show that such a verification is not possible for the other cosmologies of the  suite as the universal approach is not accurate enough to correctly capture models with massive neutrinos and dynamical dark energy.]



To date, two emulators for the halo mass function have been presented: the  <cit.> and , which is publicly available[<https://github.com/AemulusProject>] <cit.>. 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≤≤0.01 and dynamical dark energy through the w_0-w_a parametrization (see Section <ref>). 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  covers a six-dimensional cosmological parameter space. As opposed to our emulator, it is constrained to =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 <ref>). For each model, 1 Gpc/h and 2 Gpc/h simulation boxes are run evolving 2048^3 particles each; only the smaller boxes are used for the halo mass function emulator. In summary, the  framework is similar to ours except for the more restrictive range of cosmologies.

The  suite of simulations <cit.> covers seven cosmological parameters; compared to our work, they impose =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+SNIa and Planck13+BAO+SNIa data and their emulator cannot produce accurate results outside of this region of the parameter space. In contrast, the  design samples the entire parameter hypercube, as described in Section <ref>. Therefore, this hypercube can be adopted as hard priors in any cosmological analysis that uses our emulator, whereas the  emulator can only be used in analyses where strong priors from CMB, BAO, and SNIa measurements are applied. The 40  simulations have similar mass resolution as ours but smaller box sizes of 1.05 Gpc/h, thus providing a less precise calibration of the number (density) of higher-mass halos. As a cross-check, we ran the  emulator for our M000 cosmology. Note that  provides SO masses  (using the rockstar halo finder and thus not the same halo definition as we do); we convert to  assuming an NFW profile <cit.> and the concentration–mass relation by <cit.>.[We use the Colossus package: <https://bitbucket.org/bdiemer/colossus>] As shown in Figure <ref> and demonstrated in <cit.>, the prediction from  is very similar to the prediction by <cit.>: it agrees with our mass function at about 10^14M_⊙/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 <cit.>.

The discrepancy between (some of) the mass function predictions below ∼ 10^14M_⊙/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.,  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^14M_⊙/h and below will need to take this effect into account.








§ THE LIMITS OF COSMOLOGICAL UNIVERSALITY OF THE HALO MASS FUNCTION



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,cosmology). In this section, we use the  simulations to investigate the limits of the assumed cosmological universality of the mass function for .

As a first step, we establish a mass function fit following the universal approach. We choose model M000 ( with massless neutrinos) as our baseline model for which we establish the universal fit; we then extrapolate this result to all other  cosmologies and compare to the simulation results.

The universal parametrization of the halo mass function reads

    dn/dM=f(σ)ρ_m/Mdlnσ^-1/dM

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

    σ^2(M,z) = 1/2π^2∫ P(k,z) Ŵ^2(kR) k^2 dk

with the Fourier transform Ŵ of the real-space top-hat window function of radius

    R = (3M/4πρ_m)^1/3.

A common parametrization of f(σ) is

    f(σ) = B[(σ/e)^-d + σ^-f]exp(-g/σ^2)

with fit parameters d, e, f, and g <cit.> together with a normalization B set by the condition ∫ f(σ)dlnσ^-1=1 which is satisfied for

    B = 2[e^d g^-d/2Γ(d/2) + g^-f/2Γ(f/2)]^-1.

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 Eq. <ref> and matter density ρ_m that enters Eq. <ref> and <ref>. Following e.g., <cit.> we set the matter density to the density of collapsing matter:

    ρ_m = ρ_CDM+ρ_b;

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

    P(k,z) = P_CDM+b(k,z) (ρ_CDM(z)+ρ_b(z)/ρ_CDM(z)+ρ_b(z)+ρ_ν(z))^2

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 Eq. <ref> to the M000 halo catalog at z=0 using the likelihood function from Eq. <ref>. 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 <ref>, 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^14M_⊙/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.



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_180m was shown to feature a more universal behavior for  cosmologies than the mass function for  <cit.>. To extend this statement beyond , <cit.> calibrated a FOF b=0.2 mass function fit using  simulations and then used a suite of 38 simulations carried out for  cosmologies to estimate that the universal approach is still valid at the ∼10% level. Additional spot checks using a simulation with a  cosmology with massive neutrinos (νΛCDM) and a simulation with additional dynamical dark energy (νDDE) were performed by <cit.>. 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 , 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 is better understood for SO masses that are defined with respect to the critical density, such as  or M_500c, which justifies using these mass definitions despite the mass function being less universal <cit.>. 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 <ref>, 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^14M_⊙/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^15M_⊙/h.









§ SUMMARY AND OUTLOOK



We present, and make publicly available[This work uses version v0.1.0 <cit.> of the code available at <https://github.com/SebastianBocquet/MiraTitanHMFemulator>. The documentation is hosted at <https://miratitanhmfemulator.readthedocs.io/>.], an emulator for the halo mass function for redshifts z≤2 and for masses ≥10^13M_⊙/h (group and cluster-scale halos) using the newly completed  suite of 111 cosmological N-body simulations. The parameter space spans eight cosmological parameters {, , , , h, n_s, , }.
Our work expands the parameter space covered by existing mass function emulators <cit.> 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  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 Gaussian Process 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 <ref>. The accuracy is a function of mass, redshift, and location in cosmological parameter space. For example, at M=10^14M_⊙/h and z<1, the error is at the percent-level for all cosmologies; at M=10^15M_⊙/h it stays below 10%.

As pointed out in Section <ref>, 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^14M_⊙/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 non-gravitational 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 <cit.>. 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 <cit.>. As a first step, our gravity-only mass function may be modified in post-processing 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% <cit.>. It remains to be seen whether such a modified mass function prediction is accurate enough for cosmological analyses of future, high-quality datasets. 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  team for sharing some of their input mass functions for cross-checking purposes. We thank the referee for a number of thoughtful comments and suggestions that have helped improve the paper.

SB 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.
SH, KH, and EL 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 Number 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.




astropy <cit.>,
GetDist <cit.>,
matplotlib <cit.>,
numpy <cit.>,
polychord <cit.>,
pyGTC <cit.>,
(py)MultiNest <cit.>,
scipy <cit.>




We present all cosmologies of the  and of the additional test models in Table <ref>.

lcccccccc
Cosmologies in the .
Model     Ω_m h^2     Ω_b h^2
    σ_8     h     n_s     w_0     w_a     Ω_ν h^2
9lMassless-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

9lMain 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

9lAdditional 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
The models are chosen according to a tesselation-based nested design strategy <cit.>. 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.



aasjournal

