The biasing of baryons on the cluster mass function and cosmological parameter estimation
Abstract
We study the effect of baryonic processes on the halo mass function in the galaxy cluster mass range using a catalogue of 153 high resolution cosmological hydrodynamical simulations performed with the AMR code ramses. We use the results of our simulations within a simple analytical model to gauge the effects of baryon physics on the halo mass function. Neglect of AGN feedback leads to a significant boost in the cluster mass function similar to that reported by other authors. However, including AGN feedback not only gives rise to systems that are similar to observed galaxy clusters, but they also reverse the global baryonic effects on the clusters. The resulting mass function is closer to the unmodified dark matter halo mass function but still contains a mass dependent bias at the 5–10% level. These effects bias measurements of the cosmological parameters, such as and . For current cluster surveys baryonic effects are within the noise for current survey volumes, but forthcoming and planned large SZ, Xray and multiwavelength surveys will be biased at the percent level by these processes. The predictions for the halo mass function including baryonic effects need to be carefully studied with larger and improved simulations. However, simulations of full cosmological boxes with the resolution we achieve and including AGN feedback are still computationally challenging.
keywords:
black hole physics – cosmology: theory – cosmology: largescale structure of Universe – galaxies: formation – galaxies: clusters: general – methods: numerical1 Introduction
In the hierarchical structure formation scenario, the abundance of dark matter halos within a representative volume of the universe is characterized by the halo mass function. Since the statistical properties of the cosmic density field are related to the underlying cosmological model, the halo mass function contains great deal of information about the cosmological parameters of our Universe. In particular, within the standard cold dark matter cosmological scenario, dark matter halos form from initial density peaks via gravitational instability. The resulting halo mass function is directly related to the statistical properties of the primordial density field. Early studies suggested that the halo mass function can be expressed as a universal function of , the rms value of the density perturbations at a particular mass scale (Press & Schechter, 1974; Efstathiou et al., 1988; Sheth & Tormen, 1999; Sheth et al., 2001).
Cosmological Nbody simulations have been extensively used to calibrate the halo mass function and to test whether this quantity is universal and how accurately it can be determined (Audit et al., 1997; Sheth et al., 2001; Jenkins et al., 2001; Reed et al., 2003; Springel, 2005; Warren et al., 2006; Reed et al., 2007; Tinker et al., 2008; Crocce et al., 2010; Reed et al., 2013). These studies demonstrated that the mass function measured in cosmological Nbody simulations deviates from universality at the 10% level. Despite the importance of precise calibrations of the halo mass function and the ability to make accurate theoretical predictions, there is another issue that needs to be considered when studying the halo mass function: the effect of baryons on halo masses. This issue might be extremely relevant for the comparison of observationally determined halo mass functions and the prediction of cosmological simulations. In fact, present and next generation surveys such as Euclid, the Dark Energy Survey or eROSITA, are expected to measure halo masses and the halo mass function with an unprecedented percent level precision.
The effect of baryonic processes on the power spectrum and on the weak gravitational lensing shear signal has been studied in detail (White, 2004; Zhan & Knox, 2004; Jing et al., 2006; Rudd et al., 2008; Guillet et al., 2010; van Daalen et al., 2011; Semboloni et al., 2011; Reddick et al., 2013). A series of recent papers (Stanek et al., 2009; Cui et al., 2012; BalagueraAntolinez & Porciani, 2013) attempts to assess the importance of baryonic processes in shaping the properties of the halo mass function. The simulations performed by Stanek et al. (2009) and Cui et al. (2012) show a significant increase of halo counts in the galaxy cluster mass range in the cases in which cooling and star formation are considered. This is caused by the condensation of baryons at the centres of dark matter halos which leads to the contraction of the matter distribution and an enhancement of their masses. However, none of the simulations used by these authors include AGN feedback, which is expected to be an important process in shaping the properties of galaxy clusters, especially of their mass distribution (Duffy et al., 2010; Teyssier et al., 2011; Martizzi et al., 2012). In this paper, we carry out a large suite of simulations that follow the formation of galaxy clusters over a range of masses. These are fully cosmological hydrodynamical simulations performed with the ramses code. The aim of this study is to explicitly study the impact of AGN feedback on halo masses and the effects of baryonic processes on the halo mass function at redshift and the subsequent biases this introduces in recovering the underlying cosmological parameters.
The paper is organized as follows. In Section 2, we describe the numerical simulations. In Section 3, we describe the analytical formalism we adopt to compute the effect of baryons on the halo mass function. In Section 4, we discuss the measurements we perform on the simulations that will be used as an input for the analytical formalism. In Section 5, we show our predictions on the halo mass function. In Section 6, we test the impact of our predictions on the estimate of cosmological parameters. In Section 7, we summarize and discuss our results.
2 The simulations
We performed a set of 153 cosmological resimulations with the ramses code (Teyssier, 2002). We assume the standard CDM cosmological scenario with matter density parameter , cosmological constant density parameter , baryonic matter density parameter , power spectrum normalization , primordial power spectrum index and Hubble constant km/s/Mpc. The cosmological parameters are summarized in Table 1. The initial conditions for our simulations were computed using the Eisenstein & Hu (1998) transfer function and the grafic++ code developed by Doug Potter (http://sourceforge.net/projects/grafic/) which is based on the original grafic code (Bertschinger, 2001). We first ran a dark matter only simulation with particle mass M/h and box size Mpc/h. The initial level of refinement is (), but 7 more levels of refinement were carried out during the runs (maximum level ).
From this cubic cosmological volume we identified dark matter halos with the AdaptaHOP algorithm (Aubert et al., 2004), using the version implemented and tested by Tweed et al. (2009). From the list of identified halos we selected a subsample of 51 halos whose total masses are M and whose neighbouring halos do not have masses larger than within a spherical region of five times their virial radius. This choice allowed us to extract high resolution initial conditions to perform resimulations of these 51 halos. Each halo was resimulated three times: the first time considering only dark matter, the second and the third time including baryons, but using different physical prescriptions for feedback. Our results are based on these 153 simulations which constitute our cluster catalog. We label the dark matter only simulation as DMO and the hydrodynamical simulations as HYDRO.
Cosmological parameters
Type  [km sMpc]  

DMO  70.4  0.809  0.963  0.728  0.272   
HYDRO  70.4  0.809  0.963  0.728  0.272  0.045 
In the 51 DMO simulations the dark matter particle mass is M/h. In the HYDRO runs the dark matter particle mass is M/h, while the baryon resolution element has a mass of M. In all our simulations we set the maximum refinement level to which corresponds to a minimum cell size kpc/h. The grid was dynamically refined using a quasiLagrangian approach: when the dark matter or baryonic mass in a cell reaches 8 times the initial mass resolution, it is split into 8 children cells. The mass and spatial resolution of our simulations is summarized in Table 2
In the HYDRO runs, gas dynamics is modeled using a secondorder unsplit Godunov scheme (Teyssier, 2002; Teyssier et al., 2006; Fromang et al., 2006) based on the HLLC Riemann solver and the MinMod slope limiter (Toro et al., 1994). We assume a perfect gas equation of state (EOS) with polytropic index . All the HYDRO runs include subgrid models for gas cooling which account for H, He and metals and that use the Sutherland & Dopita 1993 cooling function. We directly follow star formation and supernovae feedback (”delayed cooling” scheme, Stinson et al. 2006) and metal enrichment. In 51 of the HYDRO simulations we also implement AGN feedback, using a modified version of the Booth & Schaye (2009) model in which supermassive black holes (SMBHs) are modeled as sink particles and AGN feedback is provided in form of thermal energy injected in a sphere surrounding each SMBH. We label these simulations as AGNON. The other 51 HYDRO simulations do not include AGN feedback. We label these simulations as AGNOFF.
As we showed in Martizzi et al. (2012), after implementing the physical effects of AGN feedback we were able to obtain galaxy clusters with correct stellar masses and kinematics. There are a number of unconstrained free parameters in these simulations, for example, the efficiency of the AGN feedback and star formation process. A careful study of the tuning of AGN feedback models implemented in the ramses code has been performed by Dubois et al. (2012). In our case, the tuning has been performed resimulating one of the halos in our catalog (the less massive one) several times while varying the star formation efficiency and the size of the region where the AGN feedback energy is injected. We found that the model that best reproduces the relation and the central galaxy masses has star formation efficiency and size of the AGN feedback injection region equal to twice the cell size. For the AGNOFF simulations we adopt , which is close to the lower limit of the observed star formation efficiencies. Despite this low value, Agertz et al. (2011) produced realistic ”Milky Way” candidates assuming , however we expect very different results in galaxy clusters because the AGNOFF simulations are affected by gas overcooling.
Figure 1 shows the abundance matching predictions (Moster et al., 2010) for the relation between halo mass and stellar fraction associated to the central galaxy in the clusters, compared to our results. The stellar mass of the central galaxies has been measured as in Martizzi et al. (2012). Our choice of parameters for the AGNON simulations produces results in good agreement with the data. It is worth commenting on the presence of the two outliers in the AGNON case: these two objects are highly unrelaxed clusters. In one cluster the BCG is still interacting with a close companion. In the other cluster the BCG is interacting with two companions. It is likely that taking into account the mass in the companions might improve the agreement with the abundance matching prediction. We believe that the selection criterion we adopt to identify the halos used in this study does not bias the following results. As all the unrelaxed clusters, these two outliers were excluded from the analysis concerning the mass function.
3 Effect of baryonic processes on the halo mass function
From a theoretical perspective, the easiest way to define the mass of a dark matter halo is to use spherical overdensities. We can define the average density within a given radius as:
(1) 
where is the distance from the centre of the halo and is the spherically averaged enclosed mass profile of the halo. Conventional definitions of the halo mass can be obtained by requiring to be a multiple of the critical density . For example, by requiring we define the halo mass. In this paper we are interested in studying the effect of baryonic processes related to hydrodynamics and galaxy formation on halo masses, and subsequently on the halo mass function. To do so we adopt a very simple analytical formalism, very similar to the one used by BalagueraAntolinez & Porciani (2013).
Cosmological hydrodynamical simulations have been used to show that phenomena related to galaxy formation can strongly influence the internal mass distribution of dark matter halos. Baryon condensation at the centre of halos tends to increase the concentration of the total mass distribution (Gnedin et al., 2004), whereas feedback processes tend to act in the opposite direction by decreasing the concentration through dynamical processes related to the strong winds feedback can drive (Governato et al., 2010; Duffy et al., 2010; Macciò et al., 2012; Pontzen & Governato, 2012; Martizzi et al., 2012; Teyssier et al., 2013). It has been shown that the effects of baryons are not limited to changes in the inner density profile, but to global properties such as their mass, e.g. by Stanek et al. (2009) and Cui et al. (2012) and shape, e.g. Kazantzidis et al. (2004). All of these effects can have important consequences for the halo mass function.
The analytical formalism we adopt is very simple and relies on some approximations that we outline here. We will focus on halo masses, but the formalism can be easily adapted to any definition of the halo mass. We have verified that the results we obtain are qualitatively similar in the case in which is used. We find that the adoption of is more robust because the definition of halo masses is weakly influenced by the presence of structures at the outskirts of the halo. If baryonic effects are important then if we measure , the mass of a halo in a cosmological dark matter only simulation, the mass measured in a cosmological hydrodynamical simulation will be different, . is the sum of the baryonic and nonbaryonic mass. We assume the baryon fraction within to be a function of halo mass:
(2) 
which implies that is generally different from the cosmological baryon fraction . In this paper, we adopt a very simple powerlaw scaling of the baryon fraction with mass:
(3) 
This choice is convenient since the parameters and can be easily measured from simulations as well as from observational data.
Similar formulae have been used to fit observational data, e.g. by Lin et al. (2003) and Giodini et al. (2009). This has been shown to be appropriate in the range of massive groups and clusters, however significant deviations from a powerlaw formula have been reported for very massive clusters (Leauthaud et al., 2012; Lin et al., 2012; Laganá et al., 2013). In the mass range we consider in this paper a powerlaw, like the one in Equation 2, should correctly capture the scaling of the baryon fractions with halo masses.
The effect of baryonic processes on the halo mass is modeled as a relation between and . To obtain such a relation, BalagueraAntolinez & Porciani (2013) assume that the dark matter halo mass is the same, with or without the presence of baryons. In this paper, we generalize this assumption to take account of the dynamical effect of baryons. Large concentrations of baryonic material at the centre of dark matter halos are expected to generate a contraction of the mass distribution and change the amount of dark matter present in the halos. We parametrize this effect by assuming:
(4) 
Here, we have introduced the parameter to account for the fact that the presence of baryons in the halo can induce a steepening of the density distribution within the halo and possibly boost the mass within the virial radius. In the DMO case all the mass is modeled as collisionless dark matter, but we assume that a fraction is associated to baryons. Therefore, in the DMO case we assume a dark matter mass , whereas in the HYDRO case we assume . This allows us to write Equation 4 as follows:
(5) 
We find that the parameter is quite important for the correct determination of the relation and that we cannot simply assume , i.e. we cannot assume that the dark matter mass is the same in the HYDRO and the DMO case as in BalagueraAntolinez & Porciani (2013).
In the simplistic case in which the relation between baryon fraction and mass is exact, i.e. when the intrinsic scatter is assumed to be zero, it is straightforward to compute the effect of baryonic processes on the halo mass function. Let us label the mass function affected by baryonic processes as and the standard mass function measured in dark matter only simulations as . Then, if , we have that
(6) 
Equation 6 can be easily derived by the normalization condition:
In the most general case, the relation between baryon fraction and halo mass is not a precise relation due to noise and variance, i.e. it is characterized by an average relation such as the one in Equation 2 and by a scatter about the mean. This fact means that, an value is associated to an value via a stochastic process. This process is characterized by , the conditional probability density which gives the distribution of the total cluster mass for a given . In the general case, has to be computed as
(7) 
In their discussion, BalagueraAntolinez & Porciani (2013) argue that the conditional probability function is very well approximated by a lognormal distribution and we make the same assumption here:
(8) 
where the parameters and are defined as
where is the average expected value of given , and is the scatter about . The parameters and can be easily determined from the results of our simulations. BalagueraAntolinez & Porciani (2013) argue that the lognormal approximation works very well when the logscatter is a weak function of halo mass. We find that the scatter does not show significant variations in the mass range we are considering, therefore we group all halos in a single bin to compute a unique value for . In the AGNON case we find , whereas in the AGNOFF case we find .
4 Baryon fractions and halo masses
The key elements to compute the effect of baryonic processes on the halo mass function using the formalism described in the previous section are a reliable measure of the relation and of its scatter, and a model for the conditional probability function . It is important to note that the halos considered in this paper are less massive than the ones currently used for cosmological analysis, e.g. using all sky ROSAT samples.
We measured halo masses and baryon fractions of all the halos in our cluster catalog and fitted the data using Equation 2 and Equation 5 to obtain an estimate of the average scaling relation and of its scatter. All halo masses are estimated using the spherically averaged overdensity definition of Equation 1. We then follow the formalism described in the previous section to calculate the expected halo mass functions from the three different simulations. In the following analysis we only consider relaxed clusters in order to avoid biases in the mass measurements that arise from unrelaxed systems. The relaxed clusters are identified through a morphological selection based on mock Xray emissivity maps. We identify the relaxed clusters by finding those that exhibit a unique emissivity peak within 100 kpc from their center. The selection reduces our cluster sample to 25 relaxed halos. The scatter in the properties of the relaxed sample is significantly reduced, indicating that our relaxation criterion, although minimal, does the job in selecting a uniform, well relaxed population. We show the fitting parameters in Table 3, which have been obtained using a standard algorithm. The reduced chisquared value for our fits are reported in Table 3. Such values have been obtained by weighting each term of the by the square of the inverse of 10% of the data point value. The assumed error bars on each mass measurement are supposed to mimic mass calibration errors in real surveys.
Parameters of the mass model
Type  

AGNON  0.56  
AGNOFF  0.37 
The scaling of the baryon fraction with cluster mass has been measured in observed clusters by several groups (Lin et al., 2003; Gonzalez et al., 2007; Giodini et al., 2009; Andreon, 2010; Leauthaud et al., 2012; Lin et al., 2012). The debate on which technique provides the best estimate of the baryon fraction scaling with halo mass is still ongoing (weak and strong gravitational lensing, Xray masses obtained assuming hydrostatic equilibrium, constraints from stellar kinematics, etc.). In this paper, we adopt the fits to the data provided by Lin et al. (2003) and Giodini et al. (2009), which are provided in the form of power laws similar to that in Equation 3.
Figure 2 shows a comparison of the baryon fractions measured within our virialised cluster sample at redshift with the observational results by Lin et al. (2003) and Giodini et al. (2009) (renormalized to fit our choice of cosmological parameters). As expected the AGNOFF simulations show baryon fractions in excess of the observational results, an effect related to the overcooling of gas in simulations which do not include AGN feedback. The tension between observations and simulations is partially relaxed when the AGNON case is considered. There is a marginal agreement between the AGNON results and the data by Lin et al. (2003), whereas our baryon fractions are far in excess of those measured by Giodini et al. (2009). In the AGNON case, the power law normalization we find is very similar to that found by Lin et al. (2003), however their slope is steeper. The normalization found by Giodini et al. (2009) is lower than in our AGNON case, and our slope is shallower. However, since the difference between the observed relations is large, it is not clear if our fits to the simulated data are discrepant with the real universe.
The discussion in Section 3 argues that the baryon fraction scaling with mass determines the effects of baryonic processes on the halo mass function. The data shown in Figure 2 are relevant to this. In particular, the masses of individual halos are influenced by the nature of the physical processes involved in galaxy formation (Stanek et al., 2009; Cui et al., 2012). This is indeed what we find by analysing our cluster sample, as we show in Figure 3 in which we plot the ratio as a function of . This plot demonstrates by how much the masses of the halos in the HYDRO simulations deviate from those measured for the same halos in the DMO runs. We find that in simulations which do not include AGN feedback the halo masses are boosted to values that can be more than 20% higher than the corresponding DMO value. Previous results obtained with the RAMSES code show that this effect is caused by the contraction of the total mass distribution due to the high concentration of baryons at the centre of the halo (Martizzi et al., 2012). Including AGN feedback acts in the opposite direction, preventing the formation of high concentrations of baryonic mass at the centre of the halo. The net effect is that the halos in the AGNON runs have masses that are very similar to the ones measured in the DMO case, although a significant scatter appears to be present in the scaling relation. The red and black solid lines in Figure 3 show the results of our fit to the data using the model of Equation 5.
The importance of the contraction parameter can be appreciated in Figure 4 where we plot the ratio as a function of . This plot demonstrates that the slope in the relation observed in Figure 3 is set by the variation of the baryon fractions with halo mass. The ratio of dark matter masses is approximately constant with halo mass in both the AGNON and AGNOFF cases, but it shows a large scatter. The model successfully fits the results observed in our simulations. These results have two important implications: (I) contraction of dark matter induced by baryons has a relevant effect and can be observed only in high resolution simulations; (II) the contraction effect on halo masses can be described by a simple model.
We have explicitly tested the validity of the adiabatic contraction theory (Gnedin et al., 2004) to explain the values we measure in our simulations. We adopt the adiabatic contraction model we already used in Teyssier et al. (2011) and Martizzi et al. (2012), but we assume that all the baryons are concentrated at the centre of the halo (Appendix A). The adiabatic contraction model predicts in the AGNON case and in the AGNOFF case. Both values are in good agreement with the results of our simulations, but the AGNON halos are slightly undercontracted with respect to the predictions of the adiabatic contraction model, while the AGNOFF halos are slightly overcontracted.
5 The halo mass function including baryonic effects
In the following, we assume that the halo mass function in dark matter only simulations is described by the Tinker et al. (2008) formula. The masses used in the Tinker et al. (2008) formula are defined in terms of the average matter density. To match our mass definitions with those of Tinker et al. (2008) we use the relation . We use the analytical formalism described above and the measurements performed on our resimulations to calculate the effect of baryons on the halo mass function.
The comparison between the different results is summarized in Figure 5, top panel. We only show the mass range in which the baryon fraction scaling has been measured in our simulations. Figure 5 clearly shows how the boost in halo masses observed in the AGNOFF runs influences the halo mass function by increasing the number of halos in the high mass bins. The scatter in the baryon scaling relation has been taken into account by considering the values reported at the end of Section 3 and assuming Equation 8 to be valid. The effect of the scatter of the baryon scaling relation boosts the number counts even more. As already discussed, AGN feedback acts in a way that partially suppresses the boost in halo masses observed in AGNOFF simulations. The net effect is that the prediction for the halo mass function in the AGNON case is closer to the dark matter only case. Furthermore, the effect of scatter in the baryon fraction scaling relation is also boosting the number of counts in the AGNON case. This boost is simply due to the scatter which populates higher mass bins with lower mass clusters.
There is a small caveat in this determination of the scatter calculated by comparing the same clusters simulated as dark matter only and then with baryons. Some of this scatter will arise from the numerical effects associated with different simulations of a weakly chaotic system. As an extreme example, different simulation codes that start with the same initial conditions will end up with slightly different final masses for the same object. However, we expect that most of the scatter is physical and arises from the different formation history once baryons are included. For example, cooling and starformation makes accreting halos denser and able to survive deeper into the potential and taking mass further in. We confirmed that there is a correlation between the mass increase and the stellar mass fraction of the clusters. However, to truly isolate these effects, a large number of clusters with very similar masses should be simulated.
With the aim of having a more quantitative estimate of how much the mass functions predicted for the HYDRO simulations differ from the dark matter only case, we plot the relative deviation as a function of halo mass in Figure 5, bottom panel. We also show the expected Poisson error bars on the halo mass function from a survey of volume 0.5 Gpc as a grey shaded area in Figure 5. This volume is comparable to the volume of currently available galaxy cluster surveys that have been used to obtain cosmological constraints. For instance, the MaxBCG cluster catalogue (Koester et al., 2007) covers a volume of 0.17 Gpc whereas the REFLEX II catalogue (BalagueraAntolínez et al., 2011) probes a volume of 2 Gpc. However, such a volume is quite small compared to the one that will be probed by future surveys, e.g. eROSITA, with a covered skyfraction of 0.658 and redshift limit to be 2.0, will have a survey volume larger than 400 Gpc. For Euclid the survey volume is expected to be even larger. For this reason the error bars shown in Figure 5 are much larger than the Poisson noise contribution to the expected uncertainty in the determination of the mass function from future surveys. For larger surveys the Poisson error bars will be much smaller.
The meaning of the error bars in Figure 5 is worth discussing in more detail. We are making predictions for the halo mass function, which requires halo counts in a given volume. Cosmological constraints from completed and ongoing surveys are limited by mass calibration and redshift evolution systematics (e.g. Vikhlinin et al. (2009)), therefore Poisson error bars we consider for the smaller volume represent a very idealized case. However, we are not interested in the best possible strategy for cosmological parameters extraction and in the subtleties of the estimation of the error bars, but we are interested in the effect of baryons, therefore our simplified approach should be good enough to study the problem.
Firstly, we note that the predictions of the AGNOFF model are well outside the Poisson error bars. In the case in which the scatter in the baryon fraction scaling is explicitly accounted for, the AGNOFF mass function shows a 3050% deviation from the dark matter only prediction of Tinker et al. (2008). If the AGNOFF model provided the best match to reality, this would imply serious problems for the direct comparison between measured halo mass functions and the results of dark matter only simulations: baryons would need to be explicitly and carefully accounted for when comparing measurements with simulations. However, we already know that the AGNOFF model is not a good description of observed clusters because it predicts baryon fractions and galaxy masses that are too large.
The AGNON model better matches observations, and the predictions we obtain for the halo mass function are close to the dark matter only case. However, we see that the AGNON mass function can still be distinguished from the dark matter only model assuming the Poisson error bars of a very large survey. This result is confirmed also when the scatter in the baryon fraction scaling relation is accounted for. If the results of our simulations are reliable this is a particularly important prediction, since it implies that mass functions directly measured from dark matter only simulations need to be corrected to account for the effect of baryons when comparing to observations. As we discussed above, the main effect of baryons on the halo mass function can be taken into account (e.g. by using Equations 4 to 7) once the scaling of baryon fractions with mass is known and once a simple analytical model is adopted to account for adiabatic contraction of the dark matter distribution.
6 Impact on the estimate of cosmological parameters
In the previous section we showed that including baryonic effects caused the halo mass function to differ from that found in dark matter only simulations. In this section we demonstrate that significant biases can be introduced into the estimation of cosmological parameters when baryon effects are not accounted for.
We adopt the COSMOMC (Lewis & Bridle, 2002) MCMC to infer and whilst assuming a flat universe. In our analysis, we only allowed four parameters to vary: , , and . We apply the algorithm to our mass function data using a standard Tinker et al. (2008) model to fit the data. We assume Poisson noise to compute the error bars for the resulting halo mass function. The Poisson noise is generated assuming a survey of volume 0.5 Gpc, however in full sky surveys the volumes are larger than this leading to smaller error bars and tighter constraints on the cosmological parameters. To demonstrate that, we also generate Poisson noise for a larger survey of volume 500 Gpc, comparable to eROSITA and Euclid. We stress that Poisson error bars represent an idealized case, however we do not consider this fact important for the significance of our results, as already discussed in the previous section.
First, we run the MCMC on the Tinker et al. (2008) model for the cosmology in Table 1 with the aim of obtaining likelihoods for the cosmological parameters. If the halo mass function was not influenced by baryonic effects we would get the blue contours in Figure 6, i.e. a likelihood that allows to recover the proper cosmological parameters. Our simulations show that the halo mass function is influenced by baryonic effects, so if we naively applied the Tinker et al. (2008) model to the halo mass function modified by baryonic effects we would obtain a likelihood that is biased with respect to the blue contours, because the Tinker et al. (2008) model neglects baryonic effects. This biased likelihood is represented by the red contours in Figure 6, obtained by running the MCMC on our AGNON model. Figure 6 clearly demonstrates the existence of such a bias. The top panel is for a survey of volume 0.5 Gpc, the bottom panel is for a survey of volume 500 Gpc. In the smaller volume, we find a bias in the and parameters, at the and level respectively. In the larger, volume we find a bias in the and parameters, at the and level respectively. In the bottom panel, the constraints are tighter than in the top panel (due to a larger volume), so the shift is noticeable. Because the case with baryons has a slightly higher mass function, Poisson bars are smaller and hence the constraints are also tighter. In the top panel the volume is smaller and hence the Poisson bars are larger, so the difference is not easily noticeable. Given the high degeneracy between the and parameters the overall bias is higher, up to few %. This bias will be more significant if tighter constraints are obtained from full sky surveys, as shown in the bottom panel of Figure 6, and if the constraints coming from the halo mass function are combined with those obtained from other probes (CMB, SNIa, BAO). If subpercent accuracy cosmology is the aim of the next generation observational campaign, baryonic effects clearly need to be taken into account.
7 Summary and Conclusions
We identified 51 isolated dark matter halos in a cosmological dark matter only simulation and resimulated them at higher resolution using different prescriptions for galaxy formation. The total halo masses are M M, therefore they host massive groups and clusters of galaxies. In the AGNOFF model we implement standard galaxy formation recipes, including supernovae feedback. In the AGNON model we also include AGN feedback as in Booth & Schaye (2009). We measured the masses and baryon fractions of the halos in the catalog and we used our formalism to compute how the halo mass function changes under the effect of baryons. We adopted the Tinker et al. (2008) mass function as our fiducial analytical model for the halo mass function measured in dark matter only simulations. We then used a simple analytical formalism to compute the effect of baryons on the dark matter halo mass function at redshift and we applied it to obtain predictions from cosmological hydrodynamical simulations.
Our findings can be summarized in a few points:

The AGNON model provides better fits to the observations. The masses of the central galaxies in the clusters are in agreement with the values expected from abundance matching (Moster et al., 2010; Moster et al., 2013). The scaling of baryon fractions with halo mass is in relatively good agreement with the results published by Lin et al. (2003), but in excess with respect to the values measured by Giodini et al. (2009). The AGNOFF model produces central galaxies that are too massive and baryon fractions that are too high compared to the observational results. Therefore, the AGNON model produces more realistic results with respect to the case without AGN feedback.

In the AGNOFF case halo masses can be % larger than in the dark matter only runs. This is an effect of the condensation of baryonic mass at the centre of the dark matter halos, which also boosts the concentration of the total mass distribution (Martizzi et al., 2012). AGN feedback acts in the opposite direction, preventing large condensations of baryonic mass at the centre of dark matter halos. The halo masses in the AGNON runs are very close to those measured in the dark matter only case.

Halo counts in the sampled mass range are boosted by % in the AGNOFF case. In the more realistic AGNON case we found that the halo mass function for massive groups and clusters is still not consistent with the halo mass function measured in dark matter only simulations. The halo mass function in the AGNON case is boosted by 510 %. This fact implies that mass functions measured in dark matter only simulations may not be accurate enough to be directly compared to observational results.

Neglecting baryonic effects on the halo mass function leads to biases at the % level in the estimate of cosmological parameters like and . If we assume a survey volume of 0.5 Gpc, comparable to present the one of current surveys, we find a relevant bias in the and parameters, at the and level, respectively.

Modeling of the baryonic effects on the halo mass function requires to explicitly account for the contraction of the dark matter distribution under the effect of baryons. We propose an extension of the model used by BalagueraAntolinez & Porciani (2013), which uses the baryon fraction scaling with halo mass, parameterises the response of dark matter to baryonic condensations at the halo centres and predicts the effect of baryons on halo masses.
There are a few caveats that should be discussed. First of all, we measured properties only for a limited number of halos. The bias between the pure dark matter halo mass function and that found with a treatment of hydrodynamics, starformation, supernovae and AGN feedback, may be present at the few percent level. To improve these estimates we need more statistics and a more detailed exploration of the parameter space. The best approach would be to consider a full cosmological box instead of several resimulations. However, simulations of full cosmological boxes with the resolution we achieve and including AGN feedback are still computationally challenging. This might serve as a partial justification for the adoption of the semianalytical approach we follow in this paper. Furthermore, the problem of studying the baryonic effects on the halo mass functions has also been recently studied by Velliscig et al. (2014) using SPH simulations; despite the difference of their results with ours, their conclusions is still the same: baryonic effects on the halo mass function are important. The predictions for the halo mass function including baryonic effects need to be carefully studied with larger and improved simulations, since these effects will then be relevant.
Furthermore, our analytical model to correct the mass function for baryonic effects has two limitations: firstly, the baryon fraction scaling and its scatter needs to be known; secondly, a simplified model which accounts for adiabatic contraction needs to be adopted. The former limitation could be overcome using accurate observations of the baryon fraction scaling. The latter limitation represents an intrinsic theoretical issue  the adiabatic contraction theory is only an approximation to account for the response of the mass distribution to large baryon condensations. Our results seem to confirm that such an approximation is reasonable, since the features of the halos in our numerical simulations are well captured by a simple adiabatic contraction model. We also point out that alternative techniques that take into account the effect of baryons on the halo mass function are being developed by other authors (see e.g. Cusworth et al. (2013)).
Finally, we want to stress that our results have significant implications for the correct comparison of halo counts in present and next generation SZ (e.g. Planck, South Pole Telescope, Atacama Cosmology Telescope), Xray (e.g. Chandra, XMMNewton, Suzaku, eROSITA) and multiwavelength surveys (e.g. Dark Energy Survey, Euclid): the effects of baryons influence the halo mass function at the cluster scale at a level % or less, which is larger than the level of accuracy that next generation surveys should achieve. At the galaxy cluster scale, the proper modeling of baryon physics is still likely to be an issue for what concerns the comparison between the mass function measured from observations and the one measured in simulations. For current cluster surveys baryonic effects are within the noise for current survey volumes. Recent constraints placed on cosmological parameters using cluster surveys (e.g. ROSAT (BalagueraAntolínez et al., 2011; Rapetti et al., 2013), SDSS (Zu et al., 2012; Mana et al., 2013), ACT (Hasselfield et al., 2013), SPT (Reichardt et al., 2013; Benson et al., 2013), PLANCK (Planck Collaboration et al., 2013)) are not expected to be strongly influenced by baryonic effects, but forthcoming and planned large surveys will be biased at the % level by these processes.
Acknowledgments
We thank Darren Reed for very useful discussion on the topic. We also thank the anonymous referee whose comments greatly improved the quality of this paper. All the simulations were performed on the Cray XE6 cluster Monte Rosa at CSCS, Lugano, Switzerland.
References
 Andreon (2010) Andreon S., 2010, MNRAS, 407, 263
 Aubert et al. (2004) Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376
 Audit et al. (1997) Audit E., Teyssier R., Alimi J.M., 1997, A&A, 325, 439
 BalagueraAntolinez & Porciani (2013) BalagueraAntolinez A., Porciani C., 2013, Jcap, 4, 22
 BalagueraAntolínez et al. (2011) BalagueraAntolínez A., Sánchez A. G., Böhringer H., Collins C., Guzzo L., Phleps S., 2011, MNRAS, 413, 386
 Benson et al. (2013) Benson B. A., de Haan T., Dudley J. P., Reichardt C. L., Aird K. A., Andersson K., Armstrong R., Ashby M. L. N., Bautz M., Bayliss M., Bazin G., Bleem L. E., Brodwin M., Carlstrom J. E., Chang C. L. e. a., 2013, ApJ, 763, 147
 Bertschinger (2001) Bertschinger E., 2001, The Astrophysical Journal Supplement Series, 137, 1
 Booth & Schaye (2009) Booth C. M., Schaye J., 2009, Monthly Notices of the Royal Astronomical Society, 398, 53
 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
 Cusworth et al. (2013) Cusworth S. J., Kay S. T., Battye R. A., Thomas P. A., 2013, ArXiv eprints arXiv:1309.4094
 Dubois et al. (2012) Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662
 Duffy et al. (2010) Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., Battye R. A., Booth C. M., 2010, MNRAS, 405, 2161
 Efstathiou et al. (1988) Efstathiou G., Frenk C. S., White S. D. M., Davis M., 1988, MNRAS, 235, 715
 Eisenstein & Hu (1998) Eisenstein D. J., Hu W., 1998, Astrophysical Journal v.496, 496, 605
 Fromang et al. (2006) Fromang S., Hennebelle P., Teyssier R., 2006, Astronomy and Astrophysics, 457, 371
 Giodini et al. (2009) Giodini S., Pierini D., Finoguenov A., Pratt G. W., Boehringer H., Leauthaud A., Guzzo L., Aussel H., Bolzonella M., the COSMOS Collaboration 2009, The Astrophysical Journal, 703, 982
 Gnedin et al. (2004) Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, The Astrophysical Journal, 616, 16
 Gonzalez et al. (2007) Gonzalez A. H., Zaritsky D., Zabludoff A. I., 2007, The Astrophysical Journal, 666, 147
 Governato et al. (2010) Governato F., Brook C., Mayer L., Brooks A., Rhee G., Wadsley J., Jonsson P., Willman B., Stinson G., Quinn T., Madau P., 2010, Nature, 463, 203
 Guillet et al. (2010) Guillet T., Teyssier R., Colombi S., 2010, MNRAS, 405, 525
 Hasselfield et al. (2013) Hasselfield M., Hilton M., Marriage T. A., Addison G. E., Barrientos L. F., Battaglia N., Battistelli E. S., Bond J. R., Crichton D. e. a., 2013, Jcap, 7, 8
 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
 Jing et al. (2006) Jing Y. P., Zhang P., Lin W. P., Gao L., Springel V., 2006, ApJ, 640, L119
 Kazantzidis et al. (2004) Kazantzidis S., Kravtsov A. V., Zentner A. R., Allgood B., Nagai D., Moore B., 2004, ApJ, 611, L73
 Koester et al. (2007) Koester B. P., McKay T. A., Annis J., Wechsler R. H., Evrard A., Bleem L., Becker M., Johnston D., Sheldon E., Nichol R., Miller C., Scranton R. e. a., 2007, ApJ, 660, 239
 Laganá et al. (2013) Laganá T. F., Martinet N., Durret F., Lima Neto G. B., Maughan B., Zhang Y.Y., 2013, A&A, 555, A66
 Leauthaud et al. (2012) Leauthaud A., George M. R., Behroozi P. S., Bundy K., Tinker J., Wechsler R. H., Conroy C., Finoguenov A., Tanaka M., 2012, ApJ, 746, 95
 Lewis & Bridle (2002) Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511
 Lin et al. (2003) Lin Y.T., Mohr J. J., Stanford S. A., 2003, The Astrophysical Journal, 591, 749
 Lin et al. (2012) Lin Y.T., Stanford S. A., Eisenhardt P. R. M., Vikhlinin A., Maughan B. J., Kravtsov A., 2012, ApJ, 745, L3
 Macciò et al. (2012) Macciò A. V., Stinson G., Brook C. B., Wadsley J., Couchman H. M. P., Shen S., Gibson B. K., Quinn T., 2012, ApJ, 744, L9
 Mana et al. (2013) Mana A., Giannantonio T., Weller J., Hoyle B., Hütsi G., Sartoris B., 2013, MNRAS, 434, 684
 Martizzi et al. (2012) Martizzi D., Teyssier R., Moore B., 2012, MNRAS, 420, 2859
 Martizzi et al. (2012) Martizzi D., Teyssier R., Moore B., Wentz T., 2012, MNRAS, 422, 3081
 Moster et al. (2013) Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121
 Moster et al. (2010) Moster B. P., Somerville R. S., Maulbetsch C., van den Bosch F. C., Macciò A. V., Naab T., Oser L., 2010, The Astrophysical Journal, 710, 903
 Planck Collaboration et al. (2013) Planck Collaboration Ade P. A. R., Aghanim N., ArmitageCaplan C., Arnaud M., Ashdown M., AtrioBarandela F., Aumont J., Baccigalupi C., Banday A. J. e. a., 2013, ArXiv eprints arXiv:1303.5080
 Pontzen & Governato (2012) Pontzen A., Governato F., 2012, MNRAS, 421, 3464
 Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
 Rapetti et al. (2013) Rapetti D., Blake C., Allen S. W., Mantz A., Parkinson D., Beutler F., 2013, MNRAS, 432, 973
 Reddick et al. (2013) Reddick R., Tinker J., Wechsler R., Lu Y., 2013, ArXiv eprints arXiv:1306.4686
 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
 Reichardt et al. (2013) Reichardt C. L., Stalder B., Bleem L. E., Montroy T. E., Aird K. A., Andersson K., Armstrong R., Ashby M. L. N., Bautz M., Bayliss M., Bazin G., Benson B. A., Brodwin M., Carlstrom J. E. e. a., 2013, ApJ, 763, 127
 Rudd et al. (2008) Rudd D. H., Zentner A. R., Kravtsov A. V., 2008, ApJ, 672, 19
 Semboloni et al. (2011) Semboloni E., Hoekstra H., Schaye J., van Daalen M. P., McCarthy I. G., 2011, MNRAS, 417, 2020
 Sheth et al. (2001) Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1
 Sheth & Tormen (1999) Sheth R. K., Tormen G., 1999, MNRAS, 308, 119
 Springel (2005) Springel V., 2005, MNRAS, 364, 1105
 Stanek et al. (2009) Stanek R., Rudd D., Evrard A. E., 2009, MNRAS, 394, L11
 Stinson et al. (2006) Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, Monthly Notices of the Royal Astronomical Society, 373, 1074
 Sutherland & Dopita (1993) Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253
 Teyssier (2002) Teyssier R., 2002, Astronomy and Astrophysics, 385, 337
 Teyssier et al. (2006) Teyssier R., Fromang S., Dormy E., 2006, Journal of Computational Physics, 218, 44
 Teyssier et al. (2011) Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2011, MNRAS, 414, 195
 Teyssier et al. (2013) Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068
 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
 Toro et al. (1994) Toro E. F., Spruce M., Speares W., 1994, Shock Waves, 4, 25
 Tweed et al. (2009) Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647
 van Daalen et al. (2011) van Daalen M. P., Schaye J., Booth C. M., Dalla Vecchia C., 2011, MNRAS, 415, 3649
 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, ArXiv eprints 1402.4461
 Vikhlinin et al. (2009) Vikhlinin A., Burenin R. A., Ebeling H., Forman W. R., Hornstrup A., Jones C., Kravtsov A. V., Murray S. S., Nagai D., Quintana H., Voevodkin A., 2009, ApJ, 692, 1033
 Warren et al. (2006) Warren M. S., Abazajian K., Holz D. E., Teodoro L., 2006, ApJ, 646, 881
 White (2004) White M., 2004, Astroparticle Physics, 22, 211
 Zhan & Knox (2004) Zhan H., Knox L., 2004, ApJ, 616, L75
 Zu et al. (2012) Zu Y., Weinberg D. H., Rozo E., Sheldon E. S., Tinker J. L., Becker M. R., 2012, ArXiv eprints arXiv:1207.3794
Appendix A Adiabatic Contraction Model
The simplified model we adopt is based on that used in Teyssier et al. (2011). If one defines the initial radius of each dark matter shell as , then an adiabatic contraction model is able to predict its value after contraction . In our case we adopt the transformation
(9) 
where and are the cumulative mass distributions before and after contraction, respectively. The final cumulative mass distribution can be computed as
(10) 
where is the initial mass distribution in the DMO case, is the baryonic mass distribution and is the adiabatically contracted dark matter distribution. The dark mass fraction is computed as , where is the baryon fraction. Our aim is to recover the contracted dark matter profile given and . We measure from our simulations. As a simplifying assumptions for our simple calculations, we assume that all the baryons are located at the halo centre:
(11) 
where is Dirac’s delta and is the total baryonic mass.