On the impact of baryons on the halo mass function, bias, and cluster cosmology

We define the halo mass function (HMF) as the comoving density $n(M,z)$ of halos of mass $M$ at redshift $z$, so that

is the comoving number density of such halos with mass in the range $\{M,M+{\rm d}M\}$ at the same redshift. In the above relation, the multiplicity function (MF, hereafter) $f(M,z)$ gives the fraction of total mass contained within halos of mass M, while $\rho_{\rm m}$ is the mean matter density at redshift $z=0$.

The Press-Schechter formalism (Press & Schechter, 1974) predicts the following Universal expression for the MF:

where the cosmological dependence is manifested through the peak height $\nu\equiv\nu(M,z)=\delta_{c}(z)/\sigma_{R}(M,z)$, with $\delta_{c}$ the linear density contrast for spherical collapse (Bryan & Norman, 1998). The variance of the linear density fluctuations at the scale $R$, $\sigma^{2}_{R}(z)$, is defined as

where $W(k\,R)$ is the Fourier transform of the top-hat window function of a sphere of radius $R$, and $P_{\rm lin}$ is the linear matter power spectrum. For the top-hat window function, the comoving smoothing length $R$ is related to the mass scale by $M=(4\pi/3)\,R^{3}\,\rho_{mc}.$

For the comparison with simulations, the evaluation of Eq. (3) should take into account the effects introduced by the finite simulated volume and the finite grid on which initial conditions are defined. In addition, we should also account for the specific, stochastic realisation of the initial power spectrum. We therefore compute Eq. (3) using the measured power spectrum from the sampled initial conditions. Due to the large volume of our boxes, such effects collectively affect our computation of Eq. (3) by less than $0.5$ per cent in the mass ranges defined for each box. We compared the initial conditions $P_{\rm lin}$ and the fiducial $P_{\rm lin}$. The former was calculated by integrating from $k_{\rm min}=2\pi/L_{\rm box}$ to 1/2 Nyquist frequency of the meshed grid;⁴⁴4For this test specifically we have used a grid with $1024^{3}$ mesh points and not $256^{3}$ as in the rest of the paper, see Section 4. the latter was integrated over the $k$-range from $2\pi/(10L_{\rm box})$ to 10 Nyquist frequency. Given the good agreement and in order to avoid numerical issues, such as those related to the interpolation of the power spectrum measured on the initial conditions, in the following we will use in Eq. (3) the fiducial $P_{\rm lin}$, instead of its random realization.

The MF of Eq. (2) provides a qualitative prediction for the density of objects observed in simulations. In order to have an accurate description, more sophisticated HMF modelings have been proposed and calibrated against N-body simulations. For instance, Sheth & Tormen (1999); Jenkins et al. (2001); Tinker et al. (2008, 2010); Watson et al. (2011); Despali et al. (2016); Bocquet et al. (2016) followed the approach of calibrating a parametric prescription of the HMF measured in different sets of N-body simulations, while McClintock et al. (2019); Nishimichi et al. (2019); Bocquet et al. (2020) used Gaussian process regression to build HMF emulators based on the results of large sets of N-body simulations.

In this paper, we adopt the HMF expression proposed by Tinker et al. (2008),

with the specific purpose of analysing in detail the impact of baryons. Therefore, rather than calibrating accurate fitting functions to describe the HMF from the Magneticum set of simulations, we resort to an HMF functional expression, which has been already calibrated on N-body simulations, and quantify the modification of this HMF induced by baryonic effects. As such, our HMF at $z=0$ is completely described by the $5$ parameters $\{A_{0},\,a_{0},\,b_{0},\,c_{0},\,d_{0}\}$, where $A_{0}$ is a normalization factor, $\{a_{0},\,b_{0},\,c_{0}\}$ describe the HMF at low masses and $d_{0}$ determine the location of the exponential cutoff in the high mass end. For each of these parameters, we assume a power-law dependence on the expansion factor, with exponents $\{A_{z},\,a_{z},\,b_{z},\,c_{z},\,d_{z}\}$. In total, $10$ free parameters capture the mass- and redshift-dependence of the universal (i.e. cosmology-independent) MF.

Halos are biased tracers of the underlying dark matter field. To first order, their local overdensity $\delta_{h}(\textbf{r},M,z)=n(\textbf{r},M,z)/\bar{n}(M,z)-1$ can be written as a function of the matter density contrast $\delta_{m}(\textbf{r},z)$ as

where $b$ is the halo bias and $\epsilon$ is a stochastic term, that in the following we assume to be associated to shot-noise.

From Eq. (10) it follows that, for sufficiently large scales, the halo-halo, $P_{hh}$, and halo-matter power spectrum, $P_{hm}$, are written as a function of the linear matter power spectrum, $P_{mm}$, as:

where $P_{\rm SN}$ represents a shot-noise component, which is equal to the Poisson term, $P_{\rm SN}=1/\bar{n}$, under the assumption that halos provide a discrete Poisson sampling of the underlying continuous matter density field. In fact, both negative and positive corrections to Poisson shot-noise are expected respectively for low and high mass halos (Casas-Miranda et al., 2002; Hamaus et al., 2010).

Using the Peak-Background Split (PBS) (Mo & White, 1996) it is possible to obtain the halo bias $b(M,z)$ directly from the halo mass function under the assumption of a universal MF. For the Press-Schechter MF presented in Eq. (2), the PBS prediction is:

while for the adopted HMF functional form of Eq. (4) we obtain:

In Appendix A we present the key points for the derivation of the above equation.

Although the PBS provides us a fair estimation of the bias, in Tinker et al. (2010) its performance has been observed to be not better than $10-20$ per cent. In order to obtain a better fit, we will also consider the bias fitting formula introduced in Tinker et al. (2010):

For the calibration of the HMF, we assume that the number of halos $N_{i}$ with masses $M_{\rm halo}\in[M_{i},M_{i+1})$ in the snapshot at redshift $z$, follows a Poisson distribution. The adopted likelihood $\mathcal{L}(N_{i}|\boldsymbol{\theta},z)$ is therefore given by:

where $N_{i}^{\rm sim.}$ stands for the number of halos within the above mass bin for the snapshot at redshift $z$, $N_{i}(\boldsymbol{\theta},z)$ is the theoretical expectation of halos computed by integrating Eq. (1) and multiplying it by the simulated volume assuming the HMF of Eq. (4) with parameter vector $\boldsymbol{\theta}$.

The final log-likelihood is computed by summing Eq. (16) over all redshifts, mass bins, and simulations. This amount to assume that different mass bins at fixed redshift and snapshots at different redshifts are independent of each other.

From the power spectra computed from the halo distribution with masses $M_{\rm halo}\in[M_{i},M_{i+1})$ we have selected only modes with $k$ smaller than $k_{\rm min.}=0.05\,({\rm Mpc/h})^{-1}$ in order to guarantee that the linear approximation is accurate (see Appendix B for more details).

For every mass bin, the bias is measured using the ratio of the halo-matter cross-spectrum and the matter power-spectrum: $b_{i,j}^{\rm sim.}=P_{hm}(k_{j})/P_{mm}(k_{j})$ where $i$ and $j$ indexes the mass bin and the $k$-shell.

The shell-average estimation of the power spectrum for a Gaussian field realization has a $\chi^{2}$-distribution. As the number of modes $N_{k}$ inside the shell grows rapidly with $k$, this distribution converges to a Gaussian distribution due to the central limit theorem. Thus, the distribution of the ratio is approximately described by the ratio of two Gaussian-distributed variables. In Appendix C we show that this approximation is in fact valid for our analysis.

Although computing the probability distribution function of a ratio of two random variables is straightforward, its computational cost is high if one has to compute it many times, as needed in order to find a global extrema or when running a MCMC. Fortunately, for the Gaussian case the following variable transformation leads to a normal distribution in the new variable (Geary, 1930):

where $w=y/x$, $\{\mu_{y},\sigma_{y}\}$ and $\{\mu_{x},\sigma_{x}\}$ are the mean and standard deviation of the numerator ($y$) and denominator ($x$) distributions respectively, and $\Sigma(x,y)$ is their covariance. Applying Eq. (17) to the HB likelihood estimator, we get:

with the covariances computed using linear theory:

where $\boldsymbol{\theta}$ represents the vector of parameters of Eq. (15).

The final log-likelihood is obtained summing Eq. (18) in all mass bins, selected modes, and different boxes.

Before presenting the results on the HMF and HB, we discuss in details the way in which the baryons affect such quantities. This subsection is entirely based on results from the Hydro and DMO versions of Box $2$. The numerical convergence of our results has been assessed using Box $4$ — our highest-resolution simulation. Due to its limited volume, Box $4$ allows us to assess the convergence only at at low masses and redshift $z>0.137$. In this regime, the halo statistics presented in this section is entirely consistent in both simulations, meaning that any possible high-redshift effect due to limited resolution does not propagate to low-redshift. For redshift higher than unity, Box $4$ and Box $2$ are qualitatively in agreement; however, we verified that the latter predicts halos slightly more massive than the former by $5-10$ per cent. This tension could be due to the lower resolution of Box $2$ as well as the limited volume covered by Box $4$. Thus, this section’s results are expected to be accurate at the level of $5-10$ per cent for high-redshifts.

In Figure 5 we plot the number density of halos per log-interval of the halo mass, as a function of the halo mass from our simulations. The top and bottom panels show results for masses computed at $\Delta=200m$ and $200c$ halos, respectively. The errorbars are computed using the Gaussian approximation to the Poisson distribution. Notice that the error bars are for illustration purposes only since we use a Poison likelihood on the calibration, see Section 4. Form this figure we observe small but appreciable differences between the HMF for Hydro and DMO runs in the high mass end, where Hydro runs systematically form less halos at a fixed mass.

Differences on lower masses can be better appreciated in Figure 6 where we show the relative difference of the halo number in the Hydro and in the DMO runs, relative to the best-fitting function for the DMO HMF (see below). Left and right columns are for halo mass definitions given by $\Delta=200m$ and $200c$, respectively. Different different rows are for the different redshifts considered. Dashed regions represent the $68$ per cent confidence regions for our best-fitting HMF, calculated by propagating the uncertainties on the best-fit parameters (Table 4) using the covariance matrix shown in Figures 7 and 8. For the two lowest redshifts and for masses below the minimum mass cut of Box 0 (see Table 1), we plot results from Box $2$ instead of $2b$, since only the former has been run down to $z=0.2$. Although this mass regime is a extrapolation of our fit, the performance of our fit is only slightly affected.

In order to quantify the overall quality of our HMF calibration, we present in Table 3 a two-tailed test for the reduced log-likelihood. The latter is defined by using Eq. (16), summed over all mass bins, and dividing it by the number of degrees of freedom (henceforth d.o.f.). This test has similar interpretation of the frequentist reduced-$\chi^{2}$ for Gaussian likelihoods and was done as follows: for each mass definition and simulation (Hydro and DMO), we have created $1000$ synthetic catalogs by Poisson-sampling the corresponding best fit HMF. Then, for each catalog we have computed the reduced log-likelihood at the best fit point. The $p$-value reported in Table 3 is the fraction of the catalogs that have an absolute difference between the reduced log-likelihood and the corresponding mean value of the whole set, which is smaller than the difference between our best-fit and the simulation data. The similar results for all cases indicate that our calibration performs statistically well for all of them, and the $p$-values indicate that our calibration is in agreement with the data variance, with no evidence of either under or over-fitting.

For both mass definitions, at fixed mass, halos are scarcer on Hydro than on DMO runs. The relative difference presents a characteristic mass where it is minimal and grows for both larger and smaller masses. This is the result of the interplay between the mass reduction in Hydro halos, which is less pronounced at large masses, and the exponential sensitivity of the HMF at high masses. At the low mass end and low redshifts (where critical and mean cosmic densities differ mostly), the effect of using $\Delta=200c$ for halo masses is $\sim 50$ per cent higher than for $\Delta=200m$. The enhancement of the effect when comparing critical and mean thresholds at low redshift is not surprising since the $200m$ overdensity traces larger radii where the baryonic effects are smaller.

As the calibration of the HB for $\Delta=200m$ and $200c$ involved slightly different technical aspects, we separate their presentations in Sections 5.3.1 and  5.3.2, respectively.

In Figure 11 we present our best-fit estimation for the linear halo bias, compared with the predictions of the Peak-Background Split (PBS) prescription, all computed at $z=0.293$ (our lowest redshift at which we have data from all simulations). In the top and the bottom panels we show our results for halos identified with $\Delta=200m$ and $200c$, respectively. The residual of the PBS prescription with respect to the best-fit estimation is attached on the subplot of each panel as well. The shaded area represent the $68$ per cent confidence level of the measured ratio and it has been calculated propagating the uncertainties on both HMF and linear halo bias parameters.

For the DMO results the accuracy of the PBS prediction over the $\nu$ range covered by our simulations is similar to what has been observed by Tinker et al. (2010). The quick degradation of the PBS performance at low $\nu$ is caused by our high mass cut that affects the derivative of the HMF more drastically than the HMF itself closer to the extremes. At high-$\nu$ the PBS performs with very similar accuracy on both Hydro and DMO runs. However, at intermediate $\nu$ ($0.2\lesssim\log_{10}\nu\lesssim 0.4$), the PBS performance in the Hydro is markedly worse than for the DMO. For both $200c$ and $200m$ halos, the PBS model under-predict the value of the bias by $\sim 10$ per cent for the Hydro and by few percent for the DMO halos.

The worse accuracy of the PBS prediction on the halo bias for the Hydro case is again a consequence of the variation induced by baryonic effect on halo masses. Since halos in the Hydro simulations tend to be less massive than their DMO counter-part, part of the mass inside the collapsed Lagrangian patch has not been taken into account when the mass variance is computed. In other words, baryons induce a modification of the relationship between the multiplicity function and the probability of fluctuations above a threshold on the linearly evolved density field (see Appendix A for a review of the key aspects). As the Lagrangian radius of a halo tend to be smaller in the Hydro than in the DMO, it is not surprising that also the PBS prescription underestimate the bias, as the former is a monotonic growing function of the latter.

In order to study the direct effect of baryons on clustering, we present in Figure 12 the ratio of the $P_{hh}$ computed on both Hydro and DMO runs for a matched halo catalog constructed from Box $2b$ at $z=0.25$. Different colors are different mean masses of the corresponding DMO sample. The effect of baryons on the halo clustering is fully consistent with null effect. The relative difference fluctuates around $0$ with a scatter less than $0.3$ per cent for $k<0.1$ (Mpc/$h$)^-1. This demonstrates that the baryonic effect of the mass-dependent halo bias is entirely due to the variation of halo masses, while no effect is detected on the clustering pattern of matched halo catalogues.

Galaxy cluster number counts directly constrain the halo mass function and so the cosmological parameters on which it depends.⁷⁷7See Castro et al. (2016) for a study that uses observations to constrain the halo mass function itself rather than the cosmological parameters. The number of clusters expected in a survey with sky coverage $\Omega_{\rm sky}$ within the $i$-th redshift bin $\Delta z_{i}$ centered around $z_{i}$ and the $j$-th mass bin of width $\Delta M_{j}=M_{j+1}-M_{j}$ is (see, e.g., Sartoris et al., 2016):

In the above equation $\textrm{erfc}(x_{j})$ is the complementary error function, whose argument conveys information on the halo mass through

where $M_{\rm bias}$ models a possible bias in the mass estimation (not to be confused with the bias in the halo distribution) and $\sigma_{\ln M}$ is the intrinsic scatter in the relation between true and observed mass (masses are defined according to $\Delta=200m$). Following Sartoris et al. (2016), we model the latter two quantities as

The lowest mass bin at a given redshift corresponds to $M^{ob}_{j=0}(z)=M_{\rm thr}(z)$ which defines the survey selection function, i.e. the limiting value of the observed cluster mass, $M^{ob}$ for a cluster at redshift $z$ to be included in the survey. Figure 13 reports the value of this $z$-dependent limiting mass for a Euclid-like survey, with blue and red curves corresponding to a detection threshold of signal-to-noise 5 and 3, respectively  (see Sartoris et al., 2016, Fig. 2).

Furthermore, the quantity ${\rm d}V/{\rm d}z$ appearing in Eq. (21)is the cosmology-dependent comoving volume element per unit redshift interval which is given by:

with $d_{A}$ is angular diameter distance and $H(z)$ the Hubble rate at redshift $z$.

We assume Poisson errors for the cluster counts so that we can use the Cash $C$ statistics (Cash, 1979; Holder et al., 2001):

where $L_{\rm cc}$ is the cluster count likelihood, and $N_{ij}$ and $N_{ij}^{\rm obs}$ are expected and observed counts, respectively, and $K$ is a constant. Eq. (26) is valid under the assumption that different mass- and redshift-bins are uncorrelated. Testing this hypothesis goes beyond the purpose of this analysis, as it requires using a large ensemble of mock Euclid-like surveys.

The cluster catalogs discussed in Section 6.3 can also be used to calculate the power spectrum of clusters identified according to a given selection function, $P_{\rm cl}$ (Majumdar & Mohr, 2004). It is given by (see e.g. Sartoris et al., 2010):

where $P(k,z)$ is the linear power spectrum and $b_{\rm eff}(z)$ is the linear bias weighted by the HMF:

The normalization factor $\bar{n}(z)$ is the average number density of objects included in the survey at the redshift $z$:

Note that in Eq. (27) redshift space distortions (RSD) have been neglected.

The cluster power spectrum of Eq. (27) is valid for a small redshift interval centered around $z$. Observationally, it is convenient to measure the $P_{\rm cl}$ within wide redshift intervals. The $P_{\rm cl}$ averaged over the $i$-th redshift bin $\Delta z_{i}$ centered around $z_{i}$ is then (Majumdar & Mohr, 2004):

that is, the cluster power spectrum is weighted according to the square of the number density of clusters that are included in the survey at redshift $z$.

As the so-defined $P_{\rm cl}$ probes linear scales, we can assume uncorrelated Gaussian errors so that we can build the following likelihood:

where we define again the likelihood up to a constant $K^{\prime}$ and the product runs over the redshift bins $\Delta z_{i}$ centered around $z_{i}$ and the wavenumber bins $\Delta k_{j}$ centered around $k_{j}$. We adopt constant widths for the redshift bins, $\Delta z=0.2$ (see Figure 13), and for the $k$-bin, $\Delta\log(k\,\text{Mpc}/h)=0.1$ with:

A coarser redshift bins should make correlations between adjacent bins negligible and while our choice for the value of $k_{\rm max}$ should make non-linear corrections to the power spectrum negligible (see Appendix B).

In Eq. (31) the variance is given by (Scoccimarro et al., 1999):

where $V_{k}$ is the $k$-space volume of the bin, $V_{k}=4\pi k^{2}{\rm d}k$, and $V_{s}$ is the survey volume for the redshift bin $\Delta z_{i}$, which can be computed using (25): $V_{s}=\Omega_{\rm sky}(4\pi)^{-1}\!\int_{\Delta z_{i}}\!{\rm d}z\,({\rm d}V/{\rm d}z)$. As for the forecasts of this paper, we are assuming constant (in real space) window function and power spectrum, Eq. (33) is in agreement with the optimal weighting scheme of Feldman et al. (1994). Also, Eq. (33) neglects any anisotropy in the survey volume.

Galaxy clusters sample discretely the underlying matter field and the resulting shot noise has to be accounted for. Finally, in the forecasts of Section 6.3 we model the observed power spectrum via $P_{\rm cl}^{\rm obs}=\bar{P}_{\rm cl}^{\rm fid}+1/\bar{n}^{\rm fid}$.

We consider number counts forecasts for an Euclid-like cluster survey. The Euclid telescope is scheduled for launch in 2022 and will observe approximately $\Omega_{\rm sky}=15000\deg^{2}$ of the extragalactic sky (Laureijs et al., 2011). Following Sartoris et al. (2016), the photometric cluster survey to be carried out by Euclid will have a redshift-dependent limiting mass as shown in Figure 13 for a detection signal-to-noise (S/N) threshold of 3 and 5 (see also Euclid Collaboration, 2019). The shape of the selection functions is higher at $z\sim 0.2$ than at $z\sim 0.7$, while one would in principle expect progressively more and more massive clusters to be only included at higher redshift. Sartoris et al. (2016) explains this counter-intuitive behavior as due to the competing contributions of cosmic variance and Poisson noise in the contaminating field counts. As in Sartoris et al. (2016), we assume

for the fiducial values of the four nuisance parameters of Eqs. (23)–(24). We point out that we vary these nuisance parameters in the following Fisher-Matrix forecast analysis by assuming no prior knowledge on them.

The cluster count likelihood is computed over the redshift range of Figure 13 with bins of $\Delta z=0.1$. The observed mass range extends from the lowest mass limit ($M_{\rm thr}$) to $\log M/M_{\odot}\leq 16$, with $\Delta\log_{10}M/M_{\odot}=0.2$.

As we will show in the following, neglecting baryonic effects can strongly bias cosmological inference. In the following we will show both the effect of assuming the DMO HMF and HB to analyse a population of clusters whose statistical properties are defined according to the Hydro calibration, and of correcting the DMO calibrations for baryonic effects by correcting halo masses according to what shown in Figure 1. To this end, we define the correcting function $hc(M,z)$, which is defined as the inverse of the mass variation shown in Figure 1:

In the above equation $M_{\rm DMO}$ is the mass in the DMO simulation of the halo that has mass $M_{\rm Hydro}$ in the hydro simulation. We then correct the DMO number density and bias according to: