Halo mass function: Baryon impact, fitting formulae
and implications for cluster cosmology
Abstract
We use a set of hydrodynamical (Hydro) and dark matter only (DMonly) simulations to calibrate the halo mass function (HMF). We explore the impact of baryons, propose an improved parametrization for spherical overdensity masses and identify differences between our DMonly HMF and previously published HMFs. We use the Magneticum simulations, which are well suited because of their accurate treatment of baryons, high resolution, and large cosmological volumes of up to . Baryonic effects globally decrease the masses of galaxy clusters, which, at a given mass, results in a decrease of their number density. This effect vanishes at high redshift and for high masses . We perform cosmological analyses of three idealized approximations to the cluster surveys by the South Pole Telescope (SPT), Planck, and eROSITA. We pursue two main questions: (1) What is the impact of baryons? – For the SPTlike and the Plancklike samples, the impact of baryons on cosmological results is negligible. In the eROSITAlike case, however, neglecting the baryonic impact leads to an underestimate of by about , which is comparable to the expected uncertainty from eROSITA. (2) How does our DMonly HMF compare with previous work? – For the Plancklike sample, results obtained using our DMonly HMF are shifted by with respect to results obtained using the Tinker et al. (2008) fit. This suggests that using our HMF would shift results from Planck clusters toward better agreement with CMB anisotropy measurements. Finally, we discuss biases that can be introduced through inadequate HMF parametrizations that introduce false cosmological sensitivity.
keywords:
methods: numerical – galaxies: clusters: general – cosmological parameters – cosmology: theory – dark matter – largescale structure of the Universe1 Introduction
Galaxy clusters are the largest collapsed objects in the Universe. Their distribution in mass and redshift is highly sensitive to key cosmological parameters such as the matter density , or the amount of matter fluctuations in the Universe (e.g. Henry & Arnaud, 1991; White, Efstathiou & Frenk, 1993). Furthermore, they can be used to constrain models of dark energy, the cosmic growth rate, and the neutrino sector (Wang & Steinhardt, 1998; Haiman, Mohr & Holder, 2001). Catalogues from different cluster surveys have proven to be useful cosmological probes (e.g. Vikhlinin et al., 2009; Mantz et al., 2010; Rozo et al., 2010; Benson et al., 2013; Hasselfield et al., 2013; Planck Collaboration et al., 2015b; Bocquet et al., 2015; Mantz et al., 2015).
The predicted abundance of galaxy clusters is linked to the linear matter power spectrum through the halo mass function (HMF), which was first estimated analytically (Press & Schechter, 1974). Since then, numerical body simulations have been used to calibrate fitting functions (e.g. Jenkins et al., 2001; Sheth, Mo & Tormen, 2001; White, Hernquist & Springel, 2002; Reed et al., 2003; Warren et al., 2006; Lukić et al., 2007; Reed et al., 2007; Tinker et al., 2008; Crocce et al., 2010; Bhattacharya et al., 2011; Courtin et al., 2011; Angulo et al., 2012; Watson et al., 2013). Most of the above studies focus on the friendsoffriends (FoF) halo definition (Davis et al., 1985). However, real cluster samples are typically defined in terms of spherical overdensity (SO) masses. Only very few HMFs exist for different overdensity definitions (e.g. Tinker et al., 2008; Watson et al., 2013, hereafter Tinker08, Watson13), and the former has developed into the standard reference used in most cluster cosmology analyses.
For a particular HMF parametrization to be useful in cosmological studies, it is crucially important that it correctly captures both the variation in the HMF with redshift, and the sensitivity to cosmological parameters of interest such as the matter density , the dark energy density , the dark energy equation of state parameter , and . An ideal situation would be for the HMF shape parametrization to be universal, where the variation with cosmology would be entirely captured by the cosmological sensitivity of the linear power spectrum of density fluctuations. For a FoF halo definition with linking length , or SO , the HMF has been found to be approximately universal over a wide range of redshifts and cosmologies (Jenkins et al., 2001). More recently, Bhattacharya et al. (2011) have used a set of CDM and CDM simulations to investigate the dependence of the FoF HMF on cosmological parameters. Their fit is accurate to 2% for CDM, and it describes the CDM HMF to within 10%. Similar results have also been reported in Courtin et al. (2011), although with slightly larger uncertainties.
Any HMF obtained from body dark matter only simulations potentially suffers from some bias because the baryonic component of the clusters is neglected. Recently, various authors have investigated the baryonic impact on the halo HMF using hydrodynamic simulations (e.g. Cui et al., 2012; Cui, Borgani & Murante, 2014; Cusworth et al., 2014; Martizzi et al., 2014; Velliscig et al., 2014; Vogelsberger et al., 2014; Schaller et al., 2015). Their conclusions are highly sensitive to the details of the treatment of the baryonic component. For example, models without feedback from active galactic nuclei (AGN) lead to higher cluster masses (or higher abundance at fixed mass) than dark matter only simulations. Adding AGN feedback, however, leads to HMFs that are up to 20% lower than, or about equal to the dark matter only prediction. Also, these baryonic effects are stronger for low cluster masses, and hardly affect the highmass population. These shifts in the predicted HMFs nearly correspond to the level of uncertainty from current cluster abundance measurements. Therefore, studies of the baryonic impact on the HMF are extremely important for progress in cluster cosmology.
In this work, we use haloes extracted from the hydrodynamical Magneticum simulations (Dolag et al., in prep.; see also Hirschmann et al., 2014; Saro et al., 2014; Dolag, Komatsu & Sunyaev, 2015; Teklu et al., 2015). In combination with their dark matter only counterparts, these simulations allow us to investigate several key aspects of the HMF. In particular, we focus on the impact of baryons, universality for various SO definitions, and a comparison of our dark matter only HMF with previously published fits. In Section 2, we present the Magneticum simulations and describe how the cluster samples are extracted. We discuss theoretical aspects of the HMF in Section 3, where we also present a novel approach for parametrizing the HMF for SO different from . In Section 4, we present our HMFs for different SO definitions, discuss the baryon impact and investigate differences between our dark matter only HMF and previous work. The cosmological impact of all these effects is presented in Section 5. We summarize our results in Section 6, where we also present stepbystep instructions on how to use our HMF fitting formulae.
We consider the following SO definitions: (1) mean overdensity mass , which is the mass enclosed within a sphere of radius , in which the mean enclosed matter density is equal to 200 times the mean universal matter density at the cluster’s redshift, and (2) critical overdensity masses (), which are analogous to (1) but enclosed within (), and defined with respect to the critical density . The corresponding overdensities are , and . The critical density is , where is the Hubble parameter. The mean matter density is with , and where .
Box  Size  grav. softening length (kpc)  Simulation  

DM  gas  stars  ()  ()  Hydro  DMonly  
4/uhr  ✓  ✓  –^{1}^{1}1The Hydro simulation of Box4/uhr is only run to . Its DMonly counterpart contains 835 haloes at redshift .  
3/hr  ✓  ✓  
2/hr  ✓  
2b/hr  ✓  –^{2}^{2}2Box2b/hr is only run to . At this redshift, it contains 109 578 haloes.  
1/mr  ✓  ✓  
0/mr  ✓ 
2 Simulations and Cluster selection
We will refer to our hydrodynamical simulations and HMFs as “Hydro”, and to their dark matter only counterparts as “DMonly”.
2.1 The Magneticum simulations
In this work, we use a subset of the cosmological Hydro and DMonly boxes from the Magneticum simulation set (Dolag et al., in prep.) as highlighted in Table 1. The simulations are based on the parallel cosmological TreePMSPH code PGADGET3 (Springel, 2005). We use an entropyconserving formulation of SPH (Springel & Hernquist, 2002) and a higher order kernel based on the biascorrected, sixthorder Wendland kernel (Dehnen & Aly, 2012) with 295 neighbours, which together with a lowviscosity SPH scheme allows us to properly track turbulence within galaxy clusters (Dolag et al., 2005; Donnert et al., 2013; Beck et al., 2016).
We also allow for isotropic thermal conduction with of the classical Spitzer value (Dolag et al., 2004). The simulation code includes a treatment of radiative cooling, heating from a uniform, timedependent ultraviolet background and star formation with the associated feedback processes. The latter is based on a subresolution model for the multiphase structure of the interstellar medium (Springel & Hernquist, 2003).
We compute radiative cooling rates following the same procedure as presented by Wiersma, Schaye & Smith (2009), and account for the presence of an evolving ultraviolet background (Haardt & Madau, 2001). Contributions to cooling from each element have been precomputed using the publicly available CLOUDY photoionisation code (Ferland et al., 1998) for an optically thin gas in (photo)ionisation equilibrium.
Our simulations also incorporate a detailed treatment of stellar evolution and chemical enrichment following Tornatore et al. (2007), a multiphase model for starformation (Springel & Hernquist, 2003), and feedback processes associated with supernovae driven galactic winds and AGN (Springel & Hernquist, 2003; Di Matteo et al., 2008; Fabian, 2010). Additional details about the simulation code are available elsewhere (e.g. Hirschmann et al., 2014).
Initial conditions are created at redshift applying the Zel’dovich approximation^{3}^{3}3We discuss the impact of the initial conditions in Appendix B. for a spatially flat CDM cosmology with parameter values from WMAP7 (Komatsu et al., 2011): matter density , baryon density , scalar spectral index , variance in the matter field^{4}^{4}4See Equation 2 for the definition. , and Hubble constant km s Mpc. The simulation boxes have periodic boundary conditions and are initially occupied by an equal number of gas and dark matter particles. Their relative masses reflect the global baryon fraction . To minimize numerical differences between the Hydro and the DMonly simulations, we set up the DMonly simulations with an equal number of two types of collisionless particles, whose masses are equal to the mass of the dark matter and of the gas particles, respectively, of the corresponding Hydro runs.
The Magneticum Hydro simulations have been shown to correctly reproduce AGN luminosity functions (Hirschmann et al., 2014) and cluster pressure profiles (Planck Collaboration et al., 2013; McDonald et al., 2014). In Fig. 1, we compare the baryon fraction of haloes extracted from our simulations with recent observational results on galaxy cluster and group mass scales. Simulation values at (red) are higher than the simulation values at (blue), indicating that the depletion of the intracluster medium is increasing over time. Also shown are observational results at low redshift (Laganá et al., 2011; Gonzalez et al., 2013), where we update the total and stellar masses for a subset of clusters as presented in Kravtsov, Vikhlinin & Meshscheryakov (2014), and at high redshift (Chiu et al., 2016). The low redshift observations and simulations are in good agreement, but the simulations exhibit systematically larger baryon fractions than are observed for the sample of clusters at higher redshifts.
The comparison of the observations with the simulations is complicated, given that the observations rely on measured cluster masses. As noted in Bocquet et al. (2015), the Xray hydrostatic masses tend to be lower in comparison to masses obtained using velocity dispersions or through abundance matching that includes constraints from external cosmological probes; other studies using weak lensing mass constraints come to similar conclusions (e.g. von der Linden et al., 2014; Hoekstra et al., 2015). In Fig. 1 the low redshift baryon fractions rely on Xray hydrostatic masses, and the high redshift measurements adopt the abundance matching masses referred to above. In moving from hydrostatic mass to abundance matching mass, there is an increase of 44% in the binding mass and a reduction of 27% in the baryon fraction (Chiu et al., 2016). The observed and simulated baryon fractions are in agreement if we adopt these shifts as the current scale of the systematic uncertainties.
2.2 Halo selection
Haloes are initially identified through a parallel FoF algorithm with linking length . The FoF links over dark matter particles only. We then compute SO masses (for , and ) of each halo centered at the deepest potential point with the parallel SUBFIND algorithm (Springel et al., 2001; Dolag et al., 2009). In Appendix C, we discuss an alternative approach where subhaloes are identified as individual objects.
To ensure that haloes extracted from the Hydro simulations are not affected by issues related to resolution and numerical artefacts, we apply very conservative convergence criteria. For each box, and for each overdensity , we only consider haloes that contain more than dark matter particles within . We then construct catalogues applying the lower mass limits shown in Table 1. We further apply an upper mass limit that corresponds to the lower limit of the next larger box, or to for the largest boxes (see also Fig. 2). We extract cluster catalogues at seven redshifts^{5}^{5}5We use the Hydro run of Box4/uhr for redshifts . that are roughly equally spaced in cosmic time with Gyr. This time step is chosen to be larger than the typical dynamic time of a halo, and we therefore work under the assumption that there is no correlation between the different snapshots.
3 Analysis Method
We provide the theoretical background on the HMF and introduce the fitting form we adopt. We also present the method used to perform the multidimensional fits for the HMF parameters.
3.1 The halo mass function
The comoving number density of haloes of mass is
(1) 
with the mean matter density (at redshift ), and
(2) 
which is the variance of the matter density field smoothed with the Fourier transform of the realspace tophat window function of radius . The function is commonly parametrized as
(3) 
with four parameters that need to be calibrated (e.g. Jenkins et al., 2001). Here, sets the overall normalization, and are the slope and normalization of the lowmass power law, and sets the scale of a highmass exponential cutoff.
The fit parameters in will in general depend on redshift and cosmology, which must be accounted for. In the following we note that this dependence is weak in certain special cases; this is referred to as universality of the HMF.
3.2 Halo mass function for spherical overdensity masses
Many studies of the HMF are performed using the FoF technique. For a linking length , the resulting HMF is very close to being universal (Jenkins et al., 2001). However, for observational reasons, real cluster masses are measured in terms of SO masses.
3.2.1 The halo mass function for
When using a suitable SO , the above HMF fitting function is still close to being universal (Jenkins et al., 2001). Similarly, Tinker08 use as their universal mass definition, and Watson13 argue for . These overdensity definitions are all very similar; we adopt in this work.
We follow previous work and allow for departures from universality by parametrizing a possible redshift dependence of the HMF in as a power law in :
(4) 
where the subscript 0 denotes the values at redshift , and where are additional fit parameters. Note that many authors assume the cutoff scale to be constant (, e.g. Tinker08, Watson13).
3.2.2 Toward a universal mass function for and
We also wish to calibrate the HMF for , which is a convenient mass definition within Xray studies of clusters where the emission cannot easily be traced beyond , and for , which is used for measurements of cluster galaxy velocity dispersions and for weak gravitational lensing shear profiles. It is not a priori clear that one can simply use the same form of the fitting function that is valid for , and calibrate it against simulations at other overdensities, as one might miss some redshift and cosmology dependent behaviour. Remember, for example, the very different redshift evolution of and .
Tinker08 provide the HMF for a range of different . One then uses these HMFs together with to convert from critical to mean density as a function of redshift. Their approach relies on the implicit assumption that the fitting function correctly captures the behaviour for every . Watson13 provide a correction to their HMF that depends on .
We choose a novel approach, based on the assumption that the HMF is approximately universal in . Using the HMF for we then determine the HMF for as:
(5) 
The first factor now has the same functional form as the HMF for (Equations 1–4). The second factor describes the conversion between the two different SO mass definitions. This conversion is discussed below. The crucial point is that the new functional form of the HMF established in Equation 3.2.2 has the same closetouniversal properties as the HMF for . Note that, due to uncertainties and systematic scatter in the  conversion, in Equation 3.2.2 still needs to be calibrated against the simulation data.
SO masses can be converted from one to the other assuming a cluster density profile and a massconcentration relation. We use the parametrizations by Navarro, Frenk & White (1997) and Duffy et al. (2008) and establish an analytic fit for . This fit for depends on mass, redshift, and (which is involved in the overdensity conversion). The following prescription reproduces the mass conversion at the few percent level in the range , , and :
(6) 
The parameters and are functions of and redshift:
(7) 
3.3 Finite volume correction
Throughout this work, we use cluster samples produced by simulations to understand the HMF observed in the real Universe. However, there is one subtle difference that needs to be accounted for: in contrast to the Universe, every simulation box is finite in size. Therefore, we can only capture modes in the density field that are smaller than . This means that there is an upper mass limit corresponding to the longest modes, beyond which the simulations will systematically underestimate the number of objects.
We correct for this effect following the approach of previous analyses (Lukić et al., 2007; Bhattacharya et al., 2011, Watson13). Briefly, the variance of fluctuations is corrected by subtracting the variance at scales corresponding to the box size :
(8) 
where, for simplicity, we equate the spherical volume to the cubical simulation volume . However, as we apply upper mass limits to the cluster samples (see Section 2.2), the correction has negligible impact on our analysis. In fact, for each box size, the correction would become important at masses that are about 2 orders of magnitude larger than the corresponding upper mass limit we apply. We test these finite volume effects by reproducing the results presented in Section 4 without the correction; in this case, the results do not significantly change. Nevertheless, we apply the correction to each of our boxes.
3.4 Parameter estimation
We use a Bayesian likelihood approach, which allows us to correctly capture the Poisson errors on the measured number of clusters as a function of their mass and redshift. This choice differs from using (Gaussian) statistics and jackknife errors (e.g. Tinker08), or corrections to statistics to account for the Poisson errors (e.g. Watson13).
The likelihood at each point in parameter space is calculated in the following way: We calculate the matter power spectrum using the transfer function of Eisenstein & Hu (1998, 1999), taking baryonic effects into account. This is the same prescription used to set up the initial conditions of the Magneticum simulations. We evaluate the likelihood by applying Poisson statistics in logspaced mass bins of size (Cash, 1979):
(9) 
up to an arbitrary constant, and where runs over all clusters in the sample. The second term equals the total number of expected clusters. We have checked that decreasing the bin size does not change our results.
In practice, given a set of parameters , we perform the above calculation for each redshift and for each of the simulation boxes, and sum the loglikelihoods. When fitting for the HMF in this way, we are facing a problem with moderately large dimensionality (8 free parameters); we utilize the emcee^{6}^{6}6http://dan.iel.fm/emcee/current code for efficient exploration of parameter space (ForemanMackey et al., 2013). We test our fitting procedure against several mock catalogues that contain a factor 100 times more clusters than our simulation data. In these tests we recover the input values within the statistical uncertainties and conclude that our fitting method is unbiased to a level that is much smaller than the uncertainties we report.
Parameter  

DMonly  
Hydro  
DMonly  
Hydro  
DMonly  
Hydro 
4 magneticum mass functions
In Fig. 2, we show the number density of haloes from our Hydro and DMonly simulations at seven different redshifts. The error bars show the Poisson uncertainty on the measured numbers. We also show the allowed regions of the fitting functions. The characteristic lowmass power law and the steepening at high mass are clearly visible.
The impact of baryons can be better distinguished in Fig. 3, which shows the number density of clusters from both our simulations relative to the bestfitting function for Hydro. We also show the allowed regions of our fits and two external dark matter only fitting functions (Tinker08, Watson13).
4.1 Impact of baryons
Throughout the different SO definitions, the comparison of our Hydro and DMonly simulations tells an interesting story: At , the highest redshift we consider, there is essentially no difference between our Hydro and DMonly HMFs. As structure formation continues, baryonic effects become important, and at redshift , we observe that baryon depletion is important for lowmass clusters up to . At fixed mass, this reduces the number density of clusters by . Further following the redshift evolution, the effects of baryon depletion propagate up to more massive clusters. For low redshifts , and for masses up to , our Hydro HMF is low by about compared to the DMonly case, and the difference increases when going from to to . This dependence on the SO definition is consistent with the picture of baryons being ejected from the cluster central regions, as .
In the mass range between , there is essentially no difference between our Hydro and DMonly simulations; as expected, baryonic effects are negligible on these scales. For even higher masses , and for and , the Hydro simulations contain fewer clusters than DMonly. First, we check that our simulations have converged and correctly estimate the abundance of the most massive haloes. As discussed in Appendix B, we run two additional DMonly realizations of Box1/mr, initialized at higher redshifts and , and find that they agree very well with our default simulations initialized at redshift . Then, the highmass end is much better constrained for Hydro where we also use data from Box0/mr which is not available for DMonly. We therefore expect that these differences would reduce when analyzing a larger simulation box for DMonly.
4.2 Comparison with other dark matter only fits
The comparison of our simulations and fits with Tinker08 and Watson13 can be summarized as follows:

At , the fit by Tinker08 and our DMonly fit are in very good agreement. This is an important crosscheck, since both simulations assume the same physical model and are similar in volume, and the assumed functional forms of the HMF are identical.

The fit by Watson13 tends to predict less lowmass haloes than the other HMFs considered here. At the highmass end, for , Watson13 predict less haloes than Tinker08 or our fits, but, for larger SO, Watson13 roughly agrees with Tinker08.

For the larger SO and , Tinker08 predicts up to more haloes for ; this difference steadily reduces to zero for masses .
Item (iii) above is due to a different approach in halo finding: We identify haloes and extract their masses using SUBFIND. The SO masses are constructed around the point corresponding to the minimum in the gravitational potential, neglecting substructure that eventually exceeds the overdensity when going to high values like . In contrast, Tinker08 use an SO finder that reveals these subhaloes. In Appendix C, we confirm that the presence of substructure increases the abundance of objects by about on low mass scales ; the effect drops to zero for larger masses (see also Kravtsov et al., 2004). Therefore, we expect this effect to be negligible for current cluster studies.
4.3 Halo mass function fits
In Table 2, we present the best fit parameters for the fits to the data from our DMonly and Hydro simulations. We also evaluate the per degree of freedom at the bestfitting location and the corresponding probability to exceed. The latter takes values that are very close to unity for all six HMF fits, indicating that the adopted form of the HMF is suitable for fitting the data.
The HMF parameters for can be directly compared with the literature because in this case the functional form of the HMF is identical. For DMonly, we find a value of the exponential cutoff scale that is fully consistent with Tinker08 (), but significantly smaller than the in Watson13. Note that a large value of corresponds to a low cutoff scale in mass (see Equation 3), resulting in the underprediction of massive haloes by Watson13 as shown in Fig. 3. Finally, there is no evidence for redshiftevolution of the cutoff scale in the DMonly case with , but clear evidence in the case of the Hydro HMF, where .
The HMFs for and are often used in observational studies. For the latter SO, we also show the covariance matrix for the Hydro HMF parameters in Table 3. This information should be used with the bestfitting parameters from Table 2 to capture the (statistical) uncertainties related to our HMF. For reference, the uncertainty on the HMF at is .
5 Cosmological impact
There are differences between the HMFs extracted from our Hydro and DMonly simulations, and those from the literature. When used to interpret real cluster samples, the different HMFs will ultimately lead to different cosmological results. In the following, we quantify and discuss this effect. To this end, we create simulated cluster catalogues using our bestfitting Hydro HMF, and use either the Hydro, the DMonly, or literature fits to perform cosmological analyses. Because the baryonic impact on the HMF depends on mass and redshift, we expect qualitatively different shifts when using different HMFs, depending on the properties of a specific cluster survey.
We create and analyse three sets of simulated catalogues, whose properties approximately match real samples from the South Pole Telescope (SPT; Carlstrom et al., 2011), the Planck satellite (Tauber et al., 2010) and eROSITA (Predehl et al., 2014). The selection functions we assume are shown in Fig. 4 and will be discussed in more detail. All samples are defined for SO .
We use the fit method described in Section 3, but we now fit for the cosmological parameters, and keep the HMF parameters fixed. Since no covariance matrix is available for the literature HMFs we compare to, we use our bestfit parameters without uncertainties, too, in order to make a comparison on equal footing. We further show that the (statistical) uncertainty on our Hydro HMF has negligible impact on the cosmological constraints. We restrict this analysis to the parameters and , which strongly affect the measured cluster abundance. Remember that these parameters enter the HMF calculation in Equation 1 through their impact on the matter power spectrum and the matter density . The  likelihood contours from the cluster number counts experiment exhibit a characteristic, elongated degeneracy in the  plane (see Figs 5 and 6). The parameter combination is interesting because it reflects the width of this degeneracy, i.e. the direction in  space which is best constrained using clusters. We show the constraints we recover on this parameter combination, too.
In this test, we directly use the simulated cluster masses. That is, we do not include any systematic uncertainties and measurement errors related to mass estimation as one would have to do for a real cluster sample. This also means that the uncertainties we recover only represent the statistical uncertainties related to the sample size, and cannot be compared with results from observed clusters. The aim of this analysis is to investigate and quantify offsets related to the HMF, which justifies this simplified approach. For this same reason, we do not quote the errors on the recovered parameters.
The typical uncertainties on the cosmological parameters from current cluster samples are (e.g. Planck Collaboration et al., 2014; Mantz et al., 2015). We will refer to these characteristic numbers in the following.
5.1 Cosmological analysis of an SPTlike cluster sample
The SPT sample is selected through the cluster Sunyaev–Zel’dovich effect (SZE; Sunyaev & Zel’dovich, 1972) signature, and we approximate the catalogue as massselected with , and restrict to redshifts (see Fig. 4 and cf. Bleem et al., 2015). For the SPT survey of size deg, our simulated catalogue contains systems. We consider three different input cosmologies with different values of with the same in each case. A subset of the results appears in Table 4 and Fig. 5.
The results from both our Hydro and DMonly HMFs show nearly perfect agreement, indicating that the effect of baryons on the HMF is negligible in this case. This is expected, because, as previously noted, the impact of baryons is most important for lowmass clusters at low redshifts (see Section 4.1). The SPTlike survey does not probe this mass and redshift regime.
The marginalized constraints on obtained using Tinker08 are in very good agreement with the ones from our HMFs, but slightly shifts by . Due to a mild tilt in degeneracy axes, is slightly tighter when using the Tinker08 fit.
The constraints obtained on using the Watson13 fit are tighter than the ones just discussed (see Fig. 5). However, these results seem to be biased toward . For example, the preferred value recovered for the sample with input is , and we further obtain for an input value . Their assumed form of the redshift dependence of the fit parameters (Equation 4 in this work, Equations 1315 in Watson13) involves . We suspect that this parametrization introduces an implicit and spurious preference for , which is their simulation input value. We will not consider the Watson13 fit for the rest of this work.
The statistical uncertainty of our Hydro HMF is captured by the covariance matrix in Table 3. We repeat the cosmological analysis using the Hydro HMF and its parameter covariances, and infer the additional uncertainties due to the uncertainty on the HMF using quadrature addition. We find and , and conclude that the statistical uncertainties on our HMF are completely negligible for current cluster samples.
Parameter  

Input  
SPTlike sample  
Hydro  
DMonly  
Tinker et al. (2008)  
eROSITAlike sample  
Hydro  
DMonly  
Tinker et al. (2008)  
Plancklike sample  
Input  ()  
Hydro  
DMonly  
Tinker et al. (2008) 
5.2 Cosmological analysis of a Plancklike cluster sample
The Planck cluster sample is selected using the SZE, too, and extends down to redshift . However, the satellite’s beam is larger than the SPT beam, and the survey mass limit varies significantly with redshift. We mimic the Planck selection function following the sample massredshift distribution shown in Fig. 1 in Planck Collaboration et al. (2015b). We further assume a hydrostatic bias of , which then leads to the selection function we show in Fig. 4. For this exercise, we choose our input cosmology to match the values preferred by the Planck CMB anisotropy measurement (, , km s Mpc; Planck Collaboration et al., 2015a). Assuming a sky coverage of 65%, the simulated catalogue contains clusters.
The results are shown in Table 4 and Fig. 6(a). We recover very similar constraints on for our Hydro and DMonly HMFs and Tinker08. There is some spread in the constraints on , where the difference between our DMonly and Tinker08 is . We recover identical constraints on from our Hydro and DMonly fits, but using Tinker08 leads to an shift . The shift in and the induced offset in are larger than the differences in our analysis of the SPTlike sample. As discussed in Section 4.2, our HMFs and Tinker08 differ somewhat at the highest clusters masses for redshifts . Then, as shown by the selection function in Fig. 4, the Plancksample contains more massive clusters than the SPTsample which makes the Plancksample more sensitive to differences in the highmass end of the HMF.
The analysis of the real Planck cluster sample is limited by the systematic uncertainty on the overall mass scale, often parametrized by the mass bias factor. In their latest cluster cosmology analysis, the Planck collaboration adopts measurements of the mass bias from different authors (Planck Collaboration et al., 2015b). Assuming a measurement of the mass bias from groundbased weak lensing, the recovered values of are low by roughly compared to the value recovered from the Planck CMB anisotropy measurements. However, due to the combination of systematic and statistical uncertainties, the Planck CMB preferred value in  space is still compatible with the cluster measurement at the level. Our simplified analysis suggests that using our HMF instead of Tinker08 could lead to a shift in by about , which would reduce the differences between the values preferred by Planck clusters and the CMB anisotropy.
5.3 Cosmological analysis of an eROSITAlike cluster sample
The eROSITA cluster sample will be Xray flux selected and extend from redshift . For the present test, we assume a detection limit of 50 photons in the keV band with a typical exposure time of 1.6 ks. We model this selection as a combination of a redshiftdependent mass threshold , with an additional mass cut at (see Fig. 4 and compare with fig. 2 in Pillepich, Porciani & Reiprich (2012), and also Merloni et al. (2012); Borm et al. (2014)). The eROSITA fullsky catalogue simulated in this way contains clusters.
The results of the analysis of this sample appear in Table 4 and Fig. 6(b). The recovered constraints are very tight due to the large cluster sample and the fact that we do not include mass measurement uncertainties. For this sample, the constraints on agree to within . More importantly, the recovered constraints on are in good agreement between our DMonly and Tinker08, but both are low compared to the results obtained using the Hydro HMF. This is an indication that baryonic effects are indeed important for this sample. As previously discussed, baryons have their strongest impact on the HMF at low redshifts and for low masses, which is a regime that is well probed by eROSITA. Therefore, neglecting the baryonic impact for this sample and using a dark matter only HMF would lead to a bias of . This bias is of the same order as the expected uncertainty from eROSITA (Pillepich, Porciani & Reiprich, 2012), meaning that the impact of baryons on the HMF will have to be accounted for in the cosmology analysis. For the same reasons just discussed, is slightly underestimated when using either our DMonly or the Tinker08 fit.
6 Summary
We calibrate the HMF and investigate the impact of baryons using the hydrodynamic Magneticum simulations together with dark matter only counterparts. Our simulations and the halo selection are characterized by (1) a treatment of the baryonic component and of AGN feedback that is in good agreement with several observations such as baryon fractions, AGN luminosity functions (Hirschmann et al., 2014) and cluster pressure profiles (Planck Collaboration et al., 2013; McDonald et al., 2014), (2) large cosmological volumes probed by boxes of up to , which allow us to track cluster masses up to a , and (3) a conservative halo selection with dark matter particles within , minimizing potential biases related to numerical resolution. To avoid a different sampling of the initial density fluctuations, the DMonly simulations were run using two species of dark matter with masses corresponding to those of the dark matter and baryonic particles in the Hydro simulations. We extract SO masses , , and .
The presence of baryons tends to decrease the cluster masses, which – given the shape of the HMF – leads to a decrease of the expected number of objects for a given mass (see Figs 2 and 3). The number density of haloes decreases by up to for low masses and at low redshifts . At higher masses and redshifts, our Hydro and DMonly simulations agree very well. Qualitatively similar results have been recently presented in other publications (Cui, Borgani & Murante, 2014; Cusworth et al., 2014; Velliscig et al., 2014; Vogelsberger et al., 2014; Schaller et al., 2015). Martizzi et al. (2014) find a mild increase of the HMF due to baryons.
The HMF shape varies only weakly with redshift and cosmology when masses are defined either by FoF with , or for SO . Therefore, for the HMF for (), we introduce a mapping between () and as a function of mass, redshift, and , and argue that this allows us to use the universal properties of also for masses defined by (). In practice, our HMF fits are used as follows:
Note that the same approach could be used to propagate the universal behaviour of the HMF to any overdensity .
We investigate how the differences among our Hydro, DMonly and previously published dark matter only HMFs affect cosmological results from cluster abundance measurements. To this end, we simulate idealized representations of the SPT, Planck, and eROSITA surveys, assuming simplified selection schemes as shown in Fig. 4. We assume perfect knowledge of cluster masses , and do not account for any uncertainties or systematics related to massobservable relations. Therefore, the cosmological parameter uncertainties we recover here are tighter than the actual constraints that would be obtained in a comprehensive analysis of real data. Moreover, neglecting the conversion from observable to mass likely removes some cosmological dependencies. However, this test can be used as guidance in understanding the impact of differences in the HMF.
The results of these analyses can be summarized as follows (see also Figs 5 and 6 and Table 4):

For the SPTlike and Plancklike samples, the impact of baryons is negligible, and we obtain identical cosmological results using either our Hydro or DMonly HMFs.

For the SPTlike sample, results obtained using the Tinker08 fit essentially agree with results obtained using our HMFs.

The HMF by Watson13 seems to bias results toward . This may be due to their parametrization of the redshift evolution of the HMF shape parameters using , which results in a heightened and likely artificial cosmological sensitivity.

For the Plancklike sample, using our Hydro HMF instead of Tinker08 shifts the results by . This shift corresponds to about half the observed difference between the latest Planck clusters and CMB constraints; using our HMF should therefore lead to better agreement between the two probes.

The eROSITA sample extends to lower masses than the SPT and Planck catalogues. We observe an offset in the results from Hydro and DMonly, which we identify as the impact of baryons. Neglecting this effect leads to an underestimate , which is comparable to the expected overall uncertainty.
Part of the differences between the cosmological results recovered using our HMFs and using the Tinker08 could be due to different interpretations of the connection between the observed galaxy clusters and their representation as haloes in the simulations. Different assumptions lead to different halofinding approaches: We extract SO masses around the minimum in the gravitational potential in each halo using SUBFIND. In contrast, Tinker08 employ an SO finder, and also include subhaloes in their fit. These different approaches only affect the HMF for low masses and are therefore unimportant for the cosmological study of current cluster samples.
More work, both on the theoretical and on the numerical aspects of calibrating the HMF is needed to be able to fully extract the cosmological information from nearfuture cluster samples. It is important to better understand the cosmological dependencies of the fitting functions, and to construct an analytic formula whose universality – or indeed departure from universality – is well understood. Finally, a careful comparison of cluster catalogues generated from different sets of numerical simulations would be helpful to better understand the systematic uncertainties on the HMF.
Acknowledgements
We acknowledge the support of the DFG Cluster of Excellence “Origin and Structure of the Universe” and the Transregio programme TR33 “The Dark Universe”. The calculations have partially been carried out on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP) and of the Leibniz Supercomputer Center (LRZ) under the project IDs pr58we, pr83li, and pr86re. Special thanks go to LRZ for the opportunity to run the Box0 simulation within the Extreme ScaleOut Phase on the new SuperMUC Haswell extension system. We appreciate the support from the LRZ team, especially N. Hammer, when carrying out the Box0 simulation.
References
 Angulo et al. (2012) Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
 Beck et al. (2016) Beck A. M. et al., 2016, MNRAS, 455, 2110
 Benson et al. (2013) Benson B. A. et al., 2013, ApJ, 763, 147
 Bhattacharya et al. (2011) Bhattacharya S., Heitmann K., White M., Lukić Z., Wagner C., Habib S., 2011, ApJ, 732, 122
 Bleem et al. (2015) Bleem L. E. et al., 2015, ApJS, 216, 27
 Bocquet et al. (2015) Bocquet S. et al., 2015, ApJ, 799, 214
 Borm et al. (2014) Borm K., Reiprich T. H., Mohammed I., Lovisari L., 2014, A&A, 567, A65
 Carlstrom et al. (2011) Carlstrom J. E. et al., 2011, PASP, 123, 568
 Cash (1979) Cash W., 1979, ApJ, 228, 939
 Chiu et al. (2016) Chiu I. et al., 2016, MNRAS, 455, 258
 Courtin et al. (2011) Courtin J., Rasera Y., Alimi J.M., Corasaniti P.S., Boucher V., Füzfa A., 2011, MNRAS, 410, 1911
 Crocce et al. (2010) Crocce M., Fosalba P., Castander F. J., Gaztañaga E., 2010, MNRAS, 403, 1353
 Cui et al. (2012) Cui W., Borgani S., Dolag K., Murante G., Tornatore L., 2012, MNRAS, 423, 2279
 Cui, Borgani & Murante (2014) Cui W., Borgani S., Murante G., 2014, MNRAS, 441, 1769
 Cusworth et al. (2014) Cusworth S. J., Kay S. T., Battye R. A., Thomas P. A., 2014, MNRAS, 439, 2485
 Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
 Dehnen & Aly (2012) Dehnen W., Aly H., 2012, MNRAS, 425, 1068
 Di Matteo et al. (2008) Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki D., 2008, ApJ, 676, 33
 Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
 Dolag et al. (2004) Dolag K., Jubelgas M., Springel V., Borgani S., Rasia E., 2004, ApJ, 606, L97
 Dolag, Komatsu & Sunyaev (2015) Dolag K., Komatsu E., Sunyaev R., 2015, preprint (arXiv:1509.05134)
 Dolag et al. (2005) Dolag K., Vazza F., Brunetti G., Tormen G., 2005, MNRAS, 364, 753
 Donnert et al. (2013) Donnert J., Dolag K., Brunetti G., Cassano R., 2013, MNRAS, 429, 3564
 Duffy et al. (2008) Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., 2008, MNRAS, 390, L64
 Eisenstein & Hu (1998) Eisenstein D. J., Hu W., 1998, ApJ, 496, 605
 Eisenstein & Hu (1999) Eisenstein D. J., Hu W., 1999, ApJ, 511, 5
 Fabian (2010) Fabian A. C., 2010, in IAU Symposium, Vol. 267, IAU Symposium, Peterson B. M., Somerville R. S., StorchiBergmann T., eds., pp. 341–349
 Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761
 ForemanMackey et al. (2013) ForemanMackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
 Gonzalez et al. (2013) Gonzalez A. H., Sivanandam S., Zabludoff A. I., Zaritsky D., 2013, ApJ, 778, 14
 Haardt & Madau (2001) Haardt F., Madau P., 2001, in Clusters of Galaxies and the High Redshift Universe Observed in Xrays, Neumann D. M., Tran J. T. V., eds., p. 64
 Haiman, Mohr & Holder (2001) Haiman Z., Mohr J. J., Holder G. P., 2001, ApJ, 553, 545
 Hasselfield et al. (2013) Hasselfield M. et al., 2013, J. Cosmology Astropart. Phys., 7, 8
 Henry & Arnaud (1991) Henry J. P., Arnaud K. A., 1991, ApJ, 372, 410
 Hirschmann et al. (2014) Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304
 Hoekstra et al. (2015) Hoekstra H., Herbonnet R., Muzzin A., Babul A., Mahdavi A., Viola M., Cacciato M., 2015, MNRAS, 449, 685
 Jenkins et al. (2001) Jenkins A., Frenk C. S., White S. D. M., Colberg J. M., Cole S., Evrard A. E., Couchman H. M. P., Yoshida N., 2001, MNRAS, 321, 372
 Komatsu et al. (2011) Komatsu E. et al., 2011, ApJS, 192, 18
 Kravtsov, Vikhlinin & Meshscheryakov (2014) Kravtsov A., Vikhlinin A., Meshscheryakov A., 2014, preprint (arXiv:1401.7329)
 Kravtsov et al. (2004) Kravtsov A. V., Berlind A. A., Wechsler R. H., Klypin A. A., Gottlöber S., Allgood B., Primack J. R., 2004, ApJ, 609, 35
 Laganá et al. (2011) Laganá T. F., Zhang Y.Y., Reiprich T. H., Schneider P., 2011, ApJ, 743, 13
 Lukić et al. (2007) Lukić Z., Heitmann K., Habib S., Bashinsky S., Ricker P. M., 2007, ApJ, 671, 1160
 Mantz et al. (2010) Mantz A., Allen S. W., Ebeling H., Rapetti D., DrlicaWagner A., 2010, MNRAS, 406, 1773
 Mantz et al. (2015) Mantz A. B. et al., 2015, MNRAS, 446, 2205
 Martizzi et al. (2014) Martizzi D., Mohammed I., Teyssier R., Moore B., 2014, MNRAS, 440, 2290
 McDonald et al. (2014) McDonald M. et al., 2014, ApJ, 794, 67
 Merloni et al. (2012) Merloni A. et al., 2012, preprint (arXiv:1209.3114)
 Navarro, Frenk & White (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
 Pillepich, Porciani & Reiprich (2012) Pillepich A., Porciani C., Reiprich T. H., 2012, MNRAS, 422, 44
 Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A20
 Planck Collaboration et al. (2013) Planck Collaboration et al., 2013, A&A, 550, A131
 Planck Collaboration et al. (2015a) Planck Collaboration et al., 2015a, preprint (arXiv:1502.01589)
 Planck Collaboration et al. (2015b) Planck Collaboration et al., 2015b, preprint (arXiv:1502.01597)
 Predehl et al. (2014) Predehl P. et al., 2014, in Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, Vol. 9144, Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, p. 1
 Press & Schechter (1974) Press W., Schechter P., 1974, ApJ, 187, 425
 Reed et al. (2003) Reed D., Gardner J., Quinn T., Stadel J., Fardal M., Lake G., Governato F., 2003, MNRAS, 346, 565
 Reed et al. (2007) Reed D. S., Bower R., Frenk C. S., Jenkins A., Theuns T., 2007, MNRAS, 374, 2
 Reed et al. (2013) Reed D. S., Smith R. E., Potter D., Schneider A., Stadel J., Moore B., 2013, MNRAS, 431, 1866
 Rozo et al. (2010) Rozo E. et al., 2010, ApJ, 708, 645
 Saro et al. (2014) Saro A. et al., 2014, MNRAS, 440, 2610
 Schaller et al. (2015) Schaller M. et al., 2015, MNRAS, 451, 1247
 Sheth, Mo & Tormen (2001) Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1
 Springel (2005) Springel V., 2005, MNRAS, 364, 1105
 Springel & Hernquist (2002) Springel V., Hernquist L., 2002, MNRAS, 333, 649
 Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
 Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
 Sunyaev & Zel’dovich (1972) Sunyaev R. A., Zel’dovich Y. B., 1972, Comments on Astrophysics and Space Physics, 4, 173
 Tauber et al. (2010) Tauber J. A. et al., 2010, A&A, 520, A1
 Teklu et al. (2015) Teklu A. F., Remus R.S., Dolag K., Beck A. M., Burkert A., Schmidt A. S., Schulze F., Steinborn L. K., 2015, ApJ, 812, 29
 Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
 Tornatore et al. (2007) Tornatore L., Borgani S., Dolag K., Matteucci F., 2007, MNRAS, 382, 1050
 Velliscig et al. (2014) Velliscig M., van Daalen M. P., Schaye J., McCarthy I. G., Cacciato M., Le Brun A. M. C., Dalla Vecchia C., 2014, MNRAS, 442, 2641
 Vikhlinin et al. (2009) Vikhlinin A. et al., 2009, ApJ, 692, 1060
 Vogelsberger et al. (2014) Vogelsberger M. et al., 2014, MNRAS, 444, 1518
 von der Linden et al. (2014) von der Linden A. et al., 2014, MNRAS, 443, 1973
 Wang & Steinhardt (1998) Wang L., Steinhardt P. J., 1998, ApJ, 508, 483
 Warren et al. (2006) Warren M. S., Abazajian K., Holz D. E., Teodoro L., 2006, ApJ, 646, 881
 Watson et al. (2013) Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, MNRAS, 433, 1230
 White, Hernquist & Springel (2002) White M., Hernquist L., Springel V., 2002, ApJ, 579, 16
 White, Efstathiou & Frenk (1993) White S. D. M., Efstathiou G., Frenk C. S., 1993, MNRAS, 262, 1023
 Wiersma, Schaye & Smith (2009) Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99
Appendix A Spherical overdensity
Applying the method described in Section 3.2, the HMF for is
(10) 
Assuming a cluster density profile (Navarro, Frenk & White, 1997) and a massconcentration relation (Duffy et al., 2008), we establish a massdependent fit for
(11) 
where and depend on and redshift as
(12) 
and
(13) 
This fit is accurate at the few percent level in the range , , and .
Appendix B Impact of initial redshift
The simulations used in this work are initialized at redshift using the Zel’dovich approximation. Under this assumption, running a simulation from too low an initial redshift can lead to the suppression of the formation of highmass haloes (Reed et al., 2013).
We confirm that our simulations are converged by running two additional control simulations. These correspond to Box1/mr DMonly as used in the main body of this work, except that we choose higher initial redshifts, and . In Fig. 7, we compare the abundance of haloes in the three simulations at redshifts and . All three simulations agree very well within the error bars. We conclude that our simulations initialized at redshift are converged and suitable for calibrating the HMF.
Appendix C Substructure at
As discussed in Section 4.2, the Tinker08 HMF predicts more haloes for SO and at low masses than our DMonly fit, and both HMF agree much better for . This is mainly due to different assumptions about halo definition and identification. In this section, we show that our simulations are in good agreement with the lowmass end of the Tinker08 HMF at when haloes are extracted assuming the same halo definition as applied in Tinker08. The different definitions are:

In this work, we extract spherical halo masses using the minimum in the gravitational potential within each halo as the center. In this approach, a halo with a given mass for one SO definition has exactly one counterpart in any other SO definition. Another consequence is that no subhalo that would actually exceed large overdensity thresholds like is extracted within a halo.

In Tinker08 SO masses are computed for all (potentially overlapping) haloes whose centers are separated by more than . As a consequence, for high overdensities, and within massive haloes, particles in the main halo and in its subhaloes may be counted multiple times. Furthermore, a halo identified at some SO can correspond to multiple exposed haloes at higher overdensity.
We investigate the impact of substructure at in our simulations by extracting another sample of haloes adopting the Tinker08 approach. We apply a SO finder to our DMonly simulations at redshift and extract all haloes that exceed the SO ; this sample also includes subhaloes. In Fig. 8, we compare the abundance of objects extracted in this way to the abundance of main haloes as extracted by SUBFIND and used throughout this work. Both agree for masses above , but at lower masses, subhaloes contribute to the abundance by up to for . This observation is in agreement with the discussion in Tinker08 and Kravtsov et al. (2004). Indeed, Fig. 8 suggests that, for masses below , the total abundance of haloes including subhaloes extracted from our simulations agrees with the Tinker08 HMF.
We conclude that, in the lowmass regime, the differences between the HMFs presented in this work and the HMF by Tinker08 are due to different assumptions on the halo definition. We note that these differences do not affect the shape of the HMF at higher masses which are important for the cosmological implications discussed in this work. However, this also makes it clear that the choice of (sub)halo identification applied to any observed cluster sample must be consistent with the identification method applied to the simulation data. For the cluster samples observed with the SPT and Planck, these differences are negligible, but they may indeed be important when using groups selected by, for example, eROSITA.