Constraining the mass and redshift evolution of the hydrostatic mass bias using the gas mass fraction in galaxy clusters

The gas mass fraction in galaxy clusters is a convenient probe to use in cosmological studies, as it can help derive constraints on a collection of cosmological parameters. It is however subject to various effects from the baryonic physics inside galaxy clusters, which may bias the obtained cosmological constraints. Among different aspects of the baryonic physics, in this paper we focus on the impact of the hydrostatic equilibrium assumption. We analyse the hydrostatic mass bias $B$, constraining a possible mass and redshift evolution of this quantity and its impact on the cosmological constraints. To that end we consider cluster observations of the {\it Planck}-ESZ sample and evaluate the gas mass fraction using X-ray counterpart observations. We show a degeneracy between the redshift dependence of the bias and cosmological parameters. In particular we find a $3.8 \sigma$ evidence for a redshift dependence of the bias when assuming a {\it Planck} prior on $\Omega_m$. On the other hand, assuming a constant mass bias would lead to the extreme large value of $\Omega_m>0.849$. We however show that our results are entirely dependent on the cluster sample we consider. In particular, the mass and redshift trends that we find for the lowest mass-redshift and highest mass-redshift clusters of our sample are not compatible. Nevertheless, in all the analyses we find a value for the amplitude of the bias that is consistent with $B \sim 0.8$, as expected from hydrodynamical simulations and local measurements, but still in tension with the low value of $B \sim 0.6$ derived from the combination of cosmic microwave background primary anisotropies with cluster number counts.


Introduction
Galaxy clusters are the most massive gravitationally bound systems of our Universe.As such they contain a wealth of cosmological and astrophysical information.They can be used either as powerful cosmological probes (White et al. (1993), Allen et al. (2011)) or as astrophysical objects of study to better characterise the physics of the intra-cluster medium (ICM, Sarazin (1988)) and how these structures are connected to the rest of the cosmic web.Constraints on the matter density, Ω m , or the amplitude of the matter power spectrum, σ 8 , can be inferred from several galaxy cluster observables.We could, for instance, use cluster number counts, their clustering, or the properties of their gas content to constrain cosmological parameters (see e.g.Kravtsov & Borgani (2012) or Allen et al. (2011) for reviews).Among the more recent probes, we can also cite the cluster sparsity (Balmès et al. 2013).
The baryon budget of these objects is also interesting in terms of the aspect of galaxy clusters as cosmological probes and especially the hot gas of the ICM that composes the major part of the baryonic matter inside clusters (Gouin et al. 2022).Indeed the gas mass fraction of galaxy clusters, f gas , is considered to be a good proxy for the universal baryon fraction (Borgani & Kravtsov 2011) and can be used to constrain cosmological parameters including the matter density, Ω m , the Hubble parameter, h, the dark energy density, Ω DE , or the equation of state of dark energy w (see e.g.Allen et al. (2008), Holanda et al. (2020), Mantz et al. (2022) and references therein).
The gas content inside galaxy clusters is, however, also affected by baryonic physics.Such baryonic effects need to be taken into account while performing the cosmological analysis, as they introduce systematic uncertainties in the final constraints (see e.g. the discussions from McCarthy et al. (2003a), Mc-Carthy et al. (2003b), Poole et al. (2007), Hoekstra et al. (2012), Mahdavi et al. (2013), Ruan et al. (2013), Sakr et al. (2018)).
For instance, feedback mechanisms inside clusters, such as active galactic nuclei (AGN) heating, can drive gas out of the potential wells, resulting in slightly gas-depleted clusters.This depletion of galaxy clusters' gas with respect to the universal baryon fraction is accounted for by the depletion factor Υ (Eke et al. (1998), Crain et al. (2007)).This depletion factor has been thoroughly studied in hydrodynamical simulations throughout the years, in clusters (Kravtsov et al. (2005), Planelles et al. (2013), Sembolini et al. (2016), Henden et al. (2020), and references therein) as well as in filaments (Galárraga-Espinosa et al. 2022).As a result, this parameter can be very well constrained and robustly predicted in numerical simulations.Secondly, galaxy clusters are often assumed to be in hydrostatic equilibrium (HE hereafter).Nevertheless, non-thermal processes such as turbulence, bulk motions, magnetic fields or cosmic rays (Lau et al. (2009), Vazza et al. (2009), Battaglia et al. (2012), Nelson et al. (2014), Shi et al. (2015), Biffi et al. (2016)), might cause a departure from the equilibrium condition in the ICM.
Therefore, the HE assumption leads to cluster mass estimations biased toward lower values with respect to the total cluster mass.The impact of non-thermal processes, and therefore an evaluation for the bias in the cluster mass estimation, was first considered in hydrodynamical simulations (Rasia et al. 2006).Since then, a parametrization for this mass bias has been introduced also in observations, for instance when detecting clusters in Xray or millimeter wavelengths, the latter exploiting the thermal Sunyaev-Zeldovich effect (Sunyaev & Zeldovich (1972), tSZ hereafter) (see Pratt et al. (2019) for a review).We stress that this bias affects all the observables which might assume hydrostatic equilibrium when evaluating cluster masses and, thus, the gas mass fraction.Throughout the paper we define this mass bias as B = M HE /M tot .
If the depletion factor is well constrained and understood, this is not the case for the hydrostatic mass bias as its value is still under debate.In a number of analyses including weak lensing works (Weighing the Giants (WtG) (von der Linden et al. 2014), the Canadian Cluster Comparison Project (CCCP) (Hoekstra et al. 2015), Okabe & Smith (2016), Sereno et al. (2017)), tSZ number counts (Salvati et al. 2019), X-ray observations (Eckert et al. 2019), or hydrodynamical simulations (McCarthy et al. (2017), Bennett & Sijacki (2022)), this mass bias is estimated to be around B ∼ 0.8 − 0.85.The works from Planck Collaboration et al. (2014), Planck Collaboration et al. (2016) however show that to alleviate the tension on the amplitude of the matter power spectrum, σ 8 , between local and cosmic microwave background (CMB hereafter) measurements, a much lower B is needed.Indeed, when combining CMB primary anisotropies and tSZ cluster counts, CMB measurements drive the constraining power on the cosmological parameters and, thus, on the bias, favoring a bias B ∼ 0.6 − 0.65.This result was confirmed later on by Salvati et al. (2018) andPlanck Collaboration et al. (2020).In addition, the study from Salvati et al. (2018) shows that when forcing B = 0.8 while assuming a Planck cosmology, the observed cluster number counts are way below the number counts predicted using the CMB best-fit cosmological parameters.In other words, when assuming B = 0.8 with a CMB cosmology, one predicts approximately thrice as many clusters as what is actually observed.
As the precise value of the mass bias is still an open matter and has a direct impact on the accuracy and precision of the cosmological constraints deduced from galaxy clusters, we propose a new and independent measurement of this quantity.In this paper we use the gas mass fractions of 120 galaxy clusters from the Planck-ESZ sample Planck Collaboration et al. (2011) to bring robust constraints on the value of the hydrostatic bias.More importantly, we aim at studying the potentiality of variations of the bias with mass and redshift.Such studies on mass and redshift trends of B have already been carried out in works using weak lensing (von der Linden et al. ( 2014), Hoekstra et al. (2015), Smith et al. (2016), Sereno & Ettori (2017)) or tSZ number counts (Salvati et al. 2019), sometimes giving contradictory results.In this work we measure and use gas mass fractions from a sample to get new independent constraints on B as well as on its mass and redshift evolution.This study is also aimed at investigating the role that an evolution of the bias would have on the cosmological constraints we obtain from f gas data.
After describing the theoretical modelling and our cluster sample in Section 2, we detail our methods for the data analysis in Section 3. We show our results in Section 4, first focussing on the effect of assuming a varying bias on our derived cosmological constraints, then looking at the sample dependence of our results.Finally, we discuss our results in Section 5 and draw our conclusions in Section 6.Throughout the paper, we assume a reference cosmology with H 0 = 70km.s−1 .Mpc −1 , Ω m = 0.3 and Ω Λ = 0.7.

Gas fraction sample
The gas mass fraction is defined as the ratio of the gas mass over the total mass of the cluster : Using X-ray observations, the gas mass is obtained by integrating the density profile ρ(r) inside a certain radius, r, as shown in Equation 2 below Here, the density profile ρ(r) is obtained from the electron density profile n e (r), measured from X-ray observations.We have : where µ is the mean molecular weight, m p the proton mass, n e the electron number density, and n p the proton number density with n e = 1.17n p in a fully ionised gas.Using the density profile and a temperature profile, T (r), the total hydrostatic mass, M HE , can be computed by solving the hydrostatic equilibrium equation shown below : where k B is the Boltzmann constant and G is the gravitational constant.
Similarly to the gas mass and the hydrostatic mass, the cluster gas mass fraction is evaluated within a characteristic radius, determined by the radius of the mass measurements.We define this radius, R ∆ , as the radius that encloses ∆ times the critical density of the Universe, ρ c = 3H 2 /8πG.The gas content inside galaxy clusters is affected by baryonic physics and the impact of the different astrophysical processes might depend on the considered cluster radius, R ∆ .In this work, we focus on gas fractions taken at R 500 , and from now on all the quantities we consider are taken at R 500 .We note that this radius is larger than most studies using the gas fraction of galaxy clusters as a cosmological probe, which are generally carried out at R 2500 (Allen et al. (2002), Allen et al. (2008), Allen et al. (2011), Mantz et al. (2014), Holanda et al. (2020), Mantz et al. (2022)).We briefly discuss this choice of radius in Section 5.2.2.
We computed the gas fraction for the clusters in the Planck-ESZ survey (Planck Collaboration et al. 2011).In particular, we considered 120 out of 189 clusters of the Planck-ESZ sample, for which we have follow-up X-ray observations by XMM-Newton up to R 500 (here called the ESZ sample for simplicity).Our sample therefore spans a total mass range from 2.22 × 10 14 M to 1.75×10 15 M and redshift range from 0.059 to 0.546.We started from the gas and total masses of the clusters derived in Lovisari et al. (2020b).We refer to their work for the detailed analysis of the mass evaluation.We just stress here that the total cluster masses were obtained assuming hydrostatic equilibrium, as shown in Equation 4, thus inducing the presence of the hydrostatic mass bias.From these gas masses and hydrostatic masses, we computed the gas fraction following Equation 1 to find that they are within the range [0.06, 0.20].We note that there are correlations between the gas mass and hydrostatic mass measurements.These correlations induce a corrective term when computing the error bars of f gas .We check that this correlation term (provided by Lorenzo Lovisari, private communication) introduces a negligible contribution to the total error budget.We propagate the asymmetric errors from the mass measurements to obtain the uncertainty on f gas .We also note that the Malmquist bias may affect our final results; however, the discussions in Lovisari et al. (2020b) and Andrade-Santos et al. (2021) show that this effect is negligible in the case of this sample.
We give the redshifts, gas masses, hydrostatic masses and gas fractions of the clusters from our sample in Table A.1 in the appendix.We show in Figure 1 the observed gas fraction in this sample, with respect to redshift and mass.We also compare these f gas values to the universal baryon fraction Ω b /Ω m = 0.156±0.03from Planck 2018 results (Planck Collaboration et al. 2020).

Modelling of the observed gas fraction
The hydrostatic mass is used to evaluate the total cluster mass, yet it is biased low with respect to the true total mass by a factor B = M HE /M true .The measured gas mass fraction is thus : Besides the hydrostatic mass bias, the measured gas fraction depends on a variety of instrumental, astrophysical, and cosmological effects.One way of quantifying these effects is to compare the gas fraction we obtain from gas mass and hydrostatic mass measurements to the theoretical (hydrostatic) gas fraction we expect from Equation 6 below, from Allen et al. (2008).In this equation (and in the rest of the paper), the quantity noted as X re f is the quantity X in our reference cosmology : Here, K is an instrumental calibration correction.We take K = 1 ± 0.1 from Allen et al. (2008) and discuss the soundness of this assumption in Section 5.1.1.Regarding the astrophysical contributions, Υ(M, z) is the baryon depletion factor, describing how baryons in clusters are depleted with respect to the universal baryon fraction and B(M, z) is the hydrostatic mass bias we discuss above.Finally Ω b /Ω m is the universal baryon fraction, D A is the angular diameter distance, f * is the stellar fraction, and A(z) is an angular correction which we show in Equation 7: The parameter η accounts for the slope of the f gas profiles enclosed in a spherical shell.Here we take η = 0.442 from Mantz et al. (2014).With A(z) being, however, very close to one for realistic models and within our range of parameters, the value of this parameter has a negligible impact on our results.We note that this model is valid as long we are assuming self-similarity, which we do here.Deviations from self-similarity may induce supplementary dependencies.We thus discuss our hypothesis in Section 5.1.2and Appendix B.

Methods
Our purpose in this work is to use the gas mass fraction of galaxy clusters to constrain the value of the hydrostatic mass bias and, in particular, its evolution with mass and redshift.We also want to study the role of such an evolution of the bias on the subsequent cosmological constraints derived from f gas data.We recall here that all constraints derived from gas mass fraction data are obtained by comparing the measured gas fraction to the theoretical gas fraction expected from Equation 6, which is proportional to the universal baryon fraction.Besides this proportionality, all the other constraints are deduced as corrections needed to match the observed f gas in clusters to the constant Ω b /Ω m .As discussed in Section 2, one of these correction terms accounts for the baryonic effects taking place in the ICM.These baryonic effects are accounted for by the depletion factor Υ(M, z) and the hydrostatic bias B(M, z), which we are interested in.As shown in Equation 6, we cannot constrain independently B(M, z) and Υ(M, z), as the two parameters are strongly degenerated.What we have access

Bias evolution study Sample dependence of the results Reference
Parameter Prior Prior to instead is the ratio of the two quantities, Υ(M, z)/B(M, z).In order to break this degeneracy and properly constrain the bias, strong constraints on the depletion factor and its evolution with mass and redshift are required.Obtaining such results is however out of the scope of this paper, and we retrieve these constraints from hydrodynamical simulations works.
The depletion factor is known to vary with mass, however this evolution is particularly strong for groups and low mass clusters (< 2.10 14 M ) while it becomes negligible at the high masses we consider (see the discussion in Section 3.1.1 of Eckert et al. (2019) based on results from The Three Hundred Project simulations (Cui et al. 2018)).Works from Planelles et al. (2013) and Battaglia et al. (2013) also show that the depletion factor is constant with redshift when working at R 500 .Throughout the paper, we thus assume a constant depletion factor with mass and redshift Υ 0 = 0.85 ± 0.03, based on hydrodynamical simulations from Planelles et al. (2013).

Bias evolution modelling
In order to analyze a possible mass and redshift evolution of the mass bias, we consider a power law evolution for the hydrostatic bias, with pivot masses and redshifts set at the mean values of the considered cluster sample : with B 0 as the amplitude.
We chose this model for the sake of simplicity following what is done in Salvati et al. (2019), as the exact dependence of B on mass and redshift is not known.In Section 5.1.3we discuss the role of this parametrization by comparing our results with a linear evolution of B, and find results similar to those of the power law description.The complete likelihood function thus writes, assuming no deviation from self-similarity : where with σ i as the uncertainties on the gas fraction data and f gas,T h as the gas fraction expected from Equation 6.The parameter σ f is R. Wicker , M.Douspis, L. Salvati, N. Aghanim : Constraining the mass and redshift evolution of the hydrostatic mass bias using the gas mass fraction in galaxy clusters the intrinsic scatter of the data, which we treat as a free parameter.

Free cosmology study
In the first part of our analysis we simply assess the necessity to consider an evolving bias, by looking at the impact of assuming such a variation on the subsequent cosmological constraints.To do so, we compare the posterior distributions in three different cases : In the first case we let free the baryon fraction, Ω b /Ω m , the matter density, Ω m , and the parameters accounting for the variation of the bias α, β, and B 0 .We refer to this scenario as 'VB'.In the second case we let free only the cosmological parameters and fix the bias parameters to α = 0 and β = 0, resulting in a constant bias at the value B 0 .This scenario will be noted in the rest of the paper as 'CB'.Due to a degeneracy between β and Ω m which we discuss later on, we also look at our results when leaving the set of parameters (B 0 , α, β, Ω b /Ω m ) free but assuming a prior on Ω m , in order to break this degeneracy and constrain β more accurately.We call this scenario 'VB + Ω m '.The set of parameters for which we assume flat priors in this part of the study is thus (B 0 , α, β, Ω b /Ω m , Ω m , σ f ), as we show in the first column of Table 1.This work is performed on the entirety of our cluster sample, for which our mean mass and redshift are : Throughout this whole part of the analysis, we consider a prior on the total value of the bias for a certain cluster mass and redshift, taken from the Canadian Cluster Comparison Project -Multi Epoch Nearby Cluster Survey (CCCP-MENeaCS) anal-ysis Herbonnet et al. (2020) : where z CCCP = 0.189 and M CCCP = 6.24 × 10 14 M are the mean mass and redshift for the CCCP-MENeaCS sample.

Sample dependence tests
In the second part, we looked into possible sample dependencies of our results regarding the value of the bias parameters B 0 , α, and β.To that end, we focus on different subsamples within the main sample, based on mass and redshift selections.
Matching the selection from weak lensing studies such as the Comparing Masses in Literature (CoMaLit) (Sereno & Ettori 2017) or Local Cluster Substructure Survey (LoCuSS) (Smith et al. 2016) studies, we operate a redshift cut at z = 0.2, differentiating clusters that are above or below this threshold value.This choice was also motivated by the results from Salvati et al. (2019), investigating the hydrostatic mass bias from the perspective of tSZ number counts.Their study showed that the trends in the mass bias depended on the considered redshift range, with results changing when considering only clusters with z > 0.2.We also performed a mass selection, differentiating between the clusters that are above or below the median mass of the sample M med = 5.89 × 10 14 M .In summary, the samples we consider in this study are the following : The full sample of 120 clusters, with the mean mass and redshift given previously in Equation 11.The second subsample is composed of clusters with z < 0.2 and M < 5.89 × 10 14 M .We consider them in the 'low Mz' subsample, which contains 47 clusters.The mean mass and redshift are (M lowMz , z lowMz ) = (4.26× 10 14 M , 0.126) (12) Finally our third subsample is constituted of clusters with z > 0.2 and M > 5.89 × 10 14 M .We consider them in the 'high Mz' subsample, which contains 45 clusters.The mean mass and redshift are (M highMz , z highMz ) = (8.86 × 10 14 M , 0.333) (13) We did not consider the 'low z -high M' and 'high z -low M' subsamples as by construction they do not contain enough clusters (15 and 13, respectively) to obtain meaningful results.
For illustration purposes, we show in Figure 2 the binned massredshift plane of our sample, with the delimitation of the different subsamples.Inside each bin, we computed the mean value of f gas if at least one clusters is inside the bin.
To carry out this study of the sample dependence, we compared the posterior distributions obtained when running our MCMC on the three aforementioned samples independently.We note that we kept the prior on Ω m for this part of the study, and that we add a prior on Ω b /Ω m .This choice is motivated by the presence of a degeneracy between the baryon fraction and the amplitude of the bias B 0 .This degeneracy is broken in the first part of the analysis by assuming a prior on the total value of the bias.However, when trying to compare all the bias parameters between samples (including the amplitude), we do not wish to be dependent on such a prior.The universal baryon fraction being well known and constrained, such a prior is a convenient way to break the degeneracy between Ω b /Ω m and B 0 and still obtain meaningful results on the value of all the bias parameters.The set of parameters following flat priors in this part of the study is thus (B 0 , α, β, σ f ), as we show in the second column of Table 1.
In brief, we adopted the list of priors given in Table 1 to constrain the parameters used to describe our f gas data.We fit our model in Equation 6to the f gas data with an MCMC approach using the sampler emcee (Foreman- Mackey et al. 2013).We note that the prior we consider on f * = 0.015 ± 0.005 coming from Eckert et al. (2019) has close to no effect on our final results, as this term in Equation 6 is almost negligible.

Bias evolution study
In the first part of this analysis, we intend on studying the possibility of an evolution of the hydrostatic mass bias with cluster mass and redshift.To that end, we compare the cosmological constraints obtained on the full sample when considering a constant bias to those obtained when leaving the bias free to vary.Our results are summed up in Figure 3 and in Table 2.
In the VB case, namely when leaving the set of parameters (B 0 , α, β, Ω b /Ω m , Ω m , σ f ) free with a prior on the total value of B, we show that a mass-independent bias seems to be favored.Indeed, with α = −0.056± 0.037 we remain compatible with 0 within 2σ, even though the peak is slightly lower.We also show that our derived Ω b /Ω m is fully compatible with the Planck Collaboration et al. ( 2020) value, since we obtained 0.154 +0.018 −0.026 .We cannot, however, bring such constraints on β and Ω m , as these two parameters are heavily degenerated.We zoom in on this degeneracy in Figure 4 and show that higher values of β call for higher values of Ω m .This degeneracy can be explained by the fact that both parameters entail a redshift dependence on the part of the gas fraction.Indeed, Ω m intervenes in the computation of D A (z) and H(z), which entirely drive the sensitivity of f gas to cosmology.We thus argue that this degeneracy between β and Ω m (which we show here in a simple Flat-ΛCDM model) could also be observed between β and any other cosmological parameter, as long as they appear in the computation of D A (z) or H(z).As a side note, we show a slight degeneracy between the baryon fraction Ω b /Ω m and β, with a lower Ω b /Ω m implying a higher, closer to 0 β.This degeneracy is caused by the combined effect of the degeneracy between β and Ω m , and the degeneracy between Ω b /Ω m and Ω m , which is expected (shown in Figure 3).
As the matter density, Ω m , has been strongly constrained in a number of works, we chose to assume the Planck Collaboration et al. (2020) prior shown in Table 1 on this parameter to break its degeneracy with β, in the VB + Ω m scenario.The effect of this prior is negligible on the constraints on α and Ω b /Ω m , as we obtain α = −0.057± 0.038 and 0.160 +0.016 −0.025 , fully compatible with the results obtained without this prior.On the other hand, the use of this prior allows us to constrain β.We show that the hydrostatic bias seems to show a strong redshift dependence, with β = −0.64 ± 0.18.This value of β < 0 would mean a value of B decreasing with redshift, that is to say, more and more biased masses toward higher redshifts.We also note a slight degeneracy between α and β, but this is most probably due to selection effects of the sample.As the higher redshift clusters tend to mostly have higher masses (see Figure 2), a redshift trend of the bias could then be interpreted as a slight mass trend, explaining this degeneracy.These selection effects causing the degeneracy are those we explore in the second part of the analysis when looking at the sample dependence of the results (see Sect 4.2).
The effect of assuming a constant bias (α = β = 0) in the CB case does not have a strong impact on our constraints on Ω b /Ω m , which peaks just slightly below the Planck; it does remain compatible, with Ω b /Ω m = 0.140 +0.014 −0.020 .B 0 is slightly above yet compatible with the value found in the varying bias case, with B 0 = 0.842 ± 0.040.This is simply caused by the fact that for a constant bias the amplitude B 0 is now completely determined by the total value of the bias B(z CCCP , M CCCP ).On the other hand, due to the degeneracy between β and Ω m , imposing no redshift evolution of the bias requires a very high matter density, resulting in Ω m > 0.860, fully incompatible with the Planck value.Such a high matter density is not expected in current standard cosmology and is totally aberrant.We note that these results do not completely rely on our prior on the total value of the bias.Indeed, we also performed our analysis while assuming a prior from the first results of the CCCP analysis (Hoekstra et al. 2015): with ( z , M ) = (0.246, 14.83 × 10 14 h −1 M ).Our constraints when considering this prior do not change with respect to B ( z CCCP , M CCCP ) = 0.84±0.04,except a small shift on Ω b /Ω m = 0.132 +0.019 −0.024 in the CB case.As such, we show that we need to assume an evolution of the hydrostatic mass bias, at least in redshift, to properly describe our f gas data.In the rest of our study we focus on exploring the bias evolution.We thus consider only the VB + Ω m scenario, to be able to constrain β despite its degeneracy with the matter density.We also assume a Planck Collaboration et al. ( 2020) prior on Ω b /Ω m , as the universal baryon fraction is degenerated with the total value of the bias (see Equation 6).The bias parameters (B 0 , α, β) are thus left free and we chose not to use a prior on the total value of B going forward.

Sample dependence of the results
If an evolution of the bias seems necessary to properly constrain cosmological parameters using f gas data at R 500 , the strength of this evolution might differ depending on the masses and redshifts of the clusters we consider.We thus repeat the previous study, this time focusing only on the bias parameters, using the previously defined LowMz and HighMz subsamples in addition to the full sample.We recall that in this section, we are considering the VB + Ω m scenario, with an additional prior on Ω b /Ω m .The summary of our results is given in Table 3 and in Figures 5 and  6.
Our results when considering the full sample are (as expected) the same as when adopting flat priors on the cosmological parameters.The only slight exception is B 0 , peaking slightly higher due to the absence of a prior on the total bias.We find an amplitude of B 0, f ull = 0.840 ± 0.095.The parameters accounting for the bias evolution remain unchanged, with α f ull = −0.057± 0.038 and β f ull = −0.64 ± 0.18.Thus, our results remain compatible with no mass evolution of the bias but still find a strong hint for a redshift dependence.
If all subsamples provide compatible values for the amplitude with B 0,lowMz = 0.92 +0.10 −0.11 , B 0,HighMz = 0.767 ± 0.086 and B 0, f ull = 0.840 ± 0.095, this cannot be said regarding α and β.Indeed, if the study of the full sample seems to suggest no mass evolution and a mild redshift trend of the bias, we observe the reverse behavior in the HighMz subsample.With α highMz = −0.149±0.058and β HighMz = −0.08±0.23, the preferred scenario would be the one of B constant with redshift, yet decreasing with cluster mass.On the other end of the mass-redshift plane, the results show exactly the opposite evolution, in agreement with the constraints from the full sample but aggravating the trends.With α LowMz = 0.09 ± 0.11 and β LowMz = −0.995+0.44  −0.78 , we show we are fully compatible with no mass evolution of the bias, even with the posterior of α peaking slightly above 0 contrary to the other samples.More importantly, we show a strong decreasing trend of B with redshift, as we obtain β peaking close to -1.We note that the posterior distribution for this subsample is quite wider it is than for the full sample or the HighMz subsample.This is due to the smaller mass and redshift range of this sample (see Figure 2), diminishing its constraining power with respect to the other two selections.This is, however, sufficient to highlight a redshift dependence of the bias when considering the least massive clusters of our sample at the lowest redshifts.
In Figure 6, we show what these values of the bias parameters translate to in terms of the gas fraction with respect to redshift and mass.We show these fits computed at [M min , M max ] with z free and at [z min , z max ] with M free, taking into account the uncertainties at 68% and 95% c.l.We first notice what was highlighted in the contours of Figure 5, which is that the results obtained for the full sample mainly fall in between the two extreme cases of the LowMz and HighMz subsamples.Secondly we show that the incompatibility between the two smaller subsamples is fully visible in the result of the fits.In the bottom two panels showing the f gas (M) relation, the LowMz and HighMz even seem to exhibit opposite trends, similarly to what is seen in Fig- ure 5. Finally, we note an offset in the relative positions of the curves for the two subsamples, depending on mass and redshift.This is simply due to the fact that our model describes a simultaneous evolution of the bias both with mass and redshift, which happens to be non-zero in our case.
In summary, we claim to have found strong evidence for the sample dependence of our results regarding the mass and redshift evolution of the hydrostatic mass bias.Such a sample dependence had already been noted in other works studying the evolution of the bias using tSZ cluster counts, but not when using f gas data.

Discussion
As we have shown in the previous section, considering an evolution of the bias seems to be necessary to infer sensible cosmological constraints.On the other hand, the variation of the bias that we measure is very dependent on the sample we consider, as we show different trends of the bias depending on mass and redshift selections inside our main sample.We discuss these results here, trying to take into account all the systematic effects that could appear and bias our results.We also compare our findings  to previous studies focused on the bias and its evolution, from different probes.(2021).From their work, we retrieved their total masses and compare them to the XMM-Newton masses.The result of the comparison is shown in Figure 7.We performed a fit of the point cloud in the log-log space, assuming a linear model of the form log(M XMM ) = a log(M Chandra ) + b.Similarly to the model we defined in Equation 8when studying the evolution of the bias, we set a pivot to the mean value of the Chandra masses, M Chandra = 6.32 × 10 14 M .We note that the results we present below do not change when putting the pivot at the median mass instead of the mean.To perform the fit of the point cloud we resorted to an MCMC, in order to be able to account for the asymmetric measurement errors both in X and Y axis and the intrisic scatter of the data.By doing this we obtain the following relation : log(M tot,XMM ) =(0.959 ± 0.053) log(M tot,Chandra ) (15) − (0.029 ± 0.008), with an intrinsic scatter σ i = 15.2 ± 1.3%.
We show that we are fully compatible with a massindependent mass calibration bias.We still however observe an offset, with the masses from XMM-Newton being globally lower than the Chandra masses.This result had already been shown previously (see Appendix D of Schellenberger et al. (2015)).This offset is however well accounted for by the 10% uncertainty we allow in the K prior in the majority of the sample.We do, however, point out that this comparison has been carried out only on the total masses, while our study is focused on the gas mass fraction.To be able to fully rule out possible biases from instrumental effects we would need to compare the gas mass fractions obtained from both observations rather than the total masses.Unfortunately as we do not have access to the gas masses from the Chandra observations, we can only assume that the compatibility we see for M tot stays true for f gas .This assumption is sensible, as we have As it happens, the discrepancy between Chandra and XMM-Newton mass measurements originates mainly from the calibration of the temperature profiles (Schellenberger et al. 2015).As these profiles are not used in the computation of the gas masses and as the R 500 used in both studies are compatible, is it reasonable to assume M gas,Chandra /M gas,XMM−Newton ∼ 1.We thus argue that the compatibility that we see between Chandra and XMM-Newton total mass measurements holds true for the gas fractions.
We therefore assume that if calibration biases are affecting our study, they have only minor effects on our results.

Other contributions to the evolution of f gas
In this work, we consider the evolution of the gas fraction as a probe for the evolution of the hydrostatic mass bias, as well as a cosmological probe.As discussed in Sect. 2 and shown in Eq. 6 however, several physical and instrumental effects may vary with mass and redshift and play a role in the evolution of f gas .
The depletion factor Υ is fully degenerated with the hydrostatic mass bias, both regarding its value and its evolution with mass and redshift.As such, we need a robust prior on this parameter to be able to disentangle its effects on the f gas measurements from those of the bias.Following Planelles et al. (2013), here we considered a depletion at R 500 Υ = 0.85 ± 0.03, with no evolution with mass nor redshift.This value is lower than the value of Υ = 0.94 ± 0.03 from Eckert et al. (2019), however the use of different priors on the mean value of the depletion only affects the results on B 0 , and does not change the results on the mass and redshift evolution of the bias, as we show in Appendix C. Furthermore, the latter study shows no evolution of the depletion with mass or redshift.Henden et al. (2020) found a depletion Υ ∼ 0.8 in the clusters of the Fable simulations, this time again with no clear trend with mass and redshift for clusters with M 500 > 3 × 10 14 M , compatible with the prior assumed in this work.We however stress once again that both parameters are heavily degenerated and that a shift in the value of Υ will produce a similar shift in the predicted value of B 0 .As B 0 and Ω b /Ω m are also degenerated, a shift in the value of Υ might also produce a shift in the predicted value of Ω b /Ω m .However, this has no effect on our values of α and β, nor on our constraints on other cosmological parameters.
In addition, a deviation from self-similarity due to baryonic physics may cause the M tot − M gas relation to be mass-and redshift-dependent.In Lovisari et al. (2020b), deviations from self-similarity of the ESZ cluster sample are studied.In particular, the redshift evolution of the relation is compatible with 0, in agreement with the self-similar prediction.On the other hand, the slope of the relation (i.e. its mass evolution), which we call γ departs from the self similar prediction by more than 4σ, with γ = 0.802 ± 0.049.This result however was obtained on the hydrostatic masses computed from the X-ray observations and, thus, the fitted relation is M HE − M gas .We show in Appendix B that the observed evolution of f gas with mass is, consequently, most likely due to a combination of an evolution of the bias and a deviation of the true M tot − M gas relation from self-similarity.However, we cannot disentangle the two effects and acknowledge that assuming or not self-similarity may lead to changes in the value of the mass evolution of the bias (see e.g.Eckert et al. (2016) or Truong et al. (2018)).
The calibration bias, K, may evolve with cluster mass.We show in Sect.5.1.1 that when comparing Chandra and XMM-Newton masses such an evolution should be small.For this reason, any possible impact is accounted for by our prior on K.
The stellar fraction may also vary both with mass and with redshift (see e.g.Lin et al. (2012) or Henden et al. (2020)).The mean value of this parameter being negligible, such an evolution only has a minor impact on our our results and we can consider it constant in our analysis, following the prior from Eckert et al. (2019).Thus, the only two possible remaining contributions to the evolution of f gas are the hydrostatic bias, which can evolve with mass and redshift, and the cosmology, inducing a redshift evolution on the part of the gas fraction.Following our motivated choice of priors on both bias and cosmological parameters in the two parts of the analysis, we conclude that the only effect able to induce a mass and redshift evolution of f gas is the hydrostatic mass bias.

Role of the parametrization
In the previous work Wicker et al. (2022), we investigated the evolution of B with redshift, using a linear parametrization for the evolution of the bias, instead of a power law.Here, we compare the results given by this choice of parametrization to the results we got for a power law model.In this comparison of the linear case with the power law case we focus on the redshift dependence, as we did not find strong evidence for a mass evolution of the bias in the majority of our samples.We show in Figures 8 and 9 and in Table 4 that the results are qualitatively consistent between the power law and linear descriptions.Indeed, in both cases, when considering the VB + Ω m scenario, we observe a strong sample dependence of the results.As a matter of fact, we show in Figure 8 that the LowMz clusters strongly favor a non-zero slope, namely, B 1 = −0.73+0.30 −0.66 .On the high end of the mass-redshift plane, we find results that are consistent with the power law case and that are also compatible with no redshift evolution of B, as we find B 1 = −0.095± 0.14.Our estimates of B 0 are also consistent between the power law case and the linear case for each subsample respectively, as we find B 0,lowMz = 0.912 +0.051 −0.059 , B 0,highMz = 0.763 +0.044 −0.039 and B 0, f ull = 0.839 ± 0.046.When studying the VB scenario, we still observe the strong degeneracy between the term accounting for the redshift evolution of B (here, B 1 ) and Ω m .This explains why we again obtained aberrant values of Ω m when considering the CB scenario.This compatibility between the qualitative results given by both parametrizations would hint at the fact that our results are not dominated by our choice of model for the bias evolution.We show that for low redshift and low mass clusters, the bias tends to stay constant with mass, while it increases toward higher redshift.The reverse behavior is observed in the case of high masses and high redshift clusters, with a bias that is constant with redshift but slightly increasing with the mass of the objects.This sample dependence had already been noted in Salvati et al. (2019), where the authors used tSZ number counts on the cosmology sample of the PSZ2 catalog.The authors indeed noted a non-zero trend of the bias with redshift when considering their complete sample, which disappeared when considering only the clusters above z = 0.2.We note, however, that the compatibility of our results with that study regarding the mass and redshift trends depends on the considered subsample of clusters.We ob-serve the same behavior regarding the value of the constant term of the bias B 0 , with compatible results for the high redshift clusters but incompatible ones when considering the full sample.We might argue that these differences come from different choices of priors in the two studies.Indeed, when investigating sample dependencies we did not assume any prior on the bias parameters and we considered priors on cosmology, whereas Salvati et al. (2019) assumed a prior on the total value of the bias and let their cosmological parameters remain free.However, this might not have such a strong impact, since in the first part in the analysis, we let the cosmological parameters free and we put a prior on the total value of the bias; in our results we do not see significant deviations in the value of B 0 between the two cases.These differences are actually most probably due to differences between the clusters considered in the two studies, as the clusters from the ESZ sample have globally higher mass at an equivalent redshift than the PSZ2Cosmo sample considered in their study, as we show in Figure 10.If anything, this could be yet another indication of the strong sample dependence in the results we present here.
Our results are nonetheless consistent with other weak lensing studies, which include Weighing the Giants (WtG, von der Linden et al. (2014) or CCCP (Hoekstra et al. 2015), displaying a mild decreasing trend in mass for B when considering high redshift (z WtG = 0.31, z CCCP = 0.246) and high mass (M 500,WtG = 13.58×10 14 M , M 500,CCCP = 10.38×10 14 M ) clusters.More recent works such as the X-ray study X-COP (Eckert et al. 2019) also seem to show a possible mass-dependent bias, with a decreasing trend of B. The weak-lensing studies LoCuSS (Smith et al. 2016) and CoMaLit (Sereno & Ettori 2017) both find decreasing trends of B with redshift, in agreement with this work, for clusters in the mass and redshift range of our sample (z LoCuS S = 0.22, M LoCuS S = 6.8 × 10 14 M ).
On a side note, our results regarding the amplitude, B 0 , are also compatible with Lovisari et al. (2020a), where the authors measured the ratio of hydrostatic to weak lensing mass M HE /M WL inside the Planck-ESZ sample.The authors found a ratio (M HE /M WL ) ES Z = 0.74 ± 0.06, which is in agreement with our amplitude B 0,ES Z = 0.840 ± 0.095.

Discussing the choice of a sample at R 500
In this work, we focus on gas fractions taken at R 500 , which is larger than most works using f gas as a cosmological probe, carried out at R 2500 .The choice of R 2500 is generally motivated by fraction in galaxy clusters the low scatter of the gas fraction data around those radii (see Mantz et al. (2014) and references therein), allowing for more precise cosmological constraints.The scatter f gas in data is larger at R 500 , however, this inconvenient is balanced by the stability of the depletion factor Υ at this radius.Indeed, as shown by the hydrodynamical simulations from Planelles et al. (2013), the value of Υ does not vary much with the radius when it is measured around R 500 .This reduces the possibility of a biased estimation of the depletion due to incorrect determinations of R 500 .On the other hand Υ starts to decrease in the vicinity of R 2500 .As a consequence, if the radius is not properly measured, (e.g.due to erroneous estimations of the density contrast), the estimation of the depletion will be biased.Our purpose in this study being to constrain the hydrostatic mass bias and its evolution, we need a robust prior on the depletion factor, due to the degeneracy between Υ and B. A bias in the value of the depletion due to an incorrectly measured radius would then impact our results on the mean value and evolution of B.

Discussing the tension on the bias value
Finally, we highlight that the value of B 0 = 0.840 ± 0.095 we found for the full sample is in agreement with a collection of other studies, including the aforementioned weak-lensing and Xray studies, as well as with works from hydrodynamical simulations, as shows in Figure 11.The works shown in Figure 11 are those for which the bias is known at R 500 with a central value and error bars, where the mean mass and redshift of the samples are available.We show, however, that this value is incompatible with the values B 0 = 0.62 ± 0.03 from Planck Collaboration et al. ( 2020), or B 0 0.67 from Salvati et al. (2018), which are needed to reconcile local measurements of the bias with the combination of CMB and tSZ number counts measurements.Indeed, as we show in Figure 11, the tension is alleviated only for the highest redshifts and, more particularly for the highest masses, for clusters with M 10 15 M .
Similarly the X-COP study (Eckert et al. 2019) shows that assuming a bias of B = 0.58 ± 0.04 from Planck Collaboration et al. ( 2016) results in gas fractions that are way lower than expected, as they find a median gas fraction of their sample f gas,B=0.58= 0.108 ± 0.006.This would imply that the galaxy clusters from their sample are missing about a third of their gas.A low value for the bias thus seems highly unlikely given its implications fir the gas content of galaxy clusters.
We show a similar result in this work.Indeed, from Equation 6, it is possible to compute the expected universal baryon fraction for a certain value of the bias, namely :

Conclusion
We measured the gas mass fraction of galaxy clusters at R 500 and used it to constrain a possible mass and redshift evolution of the hydrostatic mass bias.To do so, we compared the cosmological constraints we obtained when assuming a varying bias to those we obtained when assuming a constant B. We show that assuming a redshift evolution seems necessary when performing a cosmological analysis using f gas data at R 500 .Indeed we show a significant degeneracy between the redshift dependence of the bias β and the matter density Ω m .This degeneracy implies that high and close to zero β push the matter density to higher values.As a consequence, assuming a constant bias implies Ω m > 0.860, which is aberrant.Forcing Ω m to have sensible values by imposing a Planck Collaboration et al. (2020) prior in turns induces β = −0.64 ± 0.17, in a 3.8σ tension with 0. We show, however, that these results are strongly dependent on the considered sample.Indeed for the least massive clusters of our sample at the lowest redshifts, we show an important decreasing trend of B with redshift and no trend with mass, with the set (α, β) = (0.09 ± 0.11, −0.995 +0.44  −0.78 ).On the other hand, for the most massive clusters at highest redshifts, we observe no variation with redshift but we see that a decreasing evolution with mass is favored, with (α, β) = (−0.149± 0.058, −0.08 ± 0.23).When we consider the full sample, the results we obtain lie between the two extremes, largely favoring a decreasing trend of B with redshift, yet remaining compatible with no mass trend of the bias as we obtain (α, β) = (−0.057± 0.038, −0.64 ± 0.18).We recall, however, that other selection effects might affect our results as the clusters at the highest redshifts are also generally the most massive.In addition, deviations from self-similarity may lead to changes in the value of the mass evolution of the bias.
In order to identify and rule out different sources of systematic effects in our study, we looked at the possibility of being subject to instrumental calibration effects.Using mass measurements of the galaxy clusters in our sample both from Chandra and XMM-Newton, we find no evidence of a bias that would significantly change our results.In this study, we assumed a powerlaw description for the evolution of the bias.We look at the effect of this choice of parametrization by comparing our results to those we obtain when assuming a linear evolution of B with redshift.We find no major difference with the power law case, as we still observe a strong degeneracy between the redshift evolution of B and Ω m , favoring an evolution of the bias with redshift when considering a Planck prior on Ω m .Furthermore, the sample dependence we observe for the power law case holds true when assuming a linear evolution of the bias.
Despite these results, we stress that our results remain compatible with a collection of X-ray, weak lensing, and hydrodynamical simulation works regarding the mean value of the bias, given our finding of B 0 = 0.840 ± 0.095.This value remains, on the other hand, in tension with the value required from the combination of CMB observations and tSZ cluster counts to alleviate the tension on σ 8 .
Finally, we recall that this work is focused on gas fractions taken at R 500 , with a goal to obtain constraints on two cosmological parameters, the universal baryon fraction, Ω b /Ω m , and the matter density of the Universe, Ω m , which are the two parameters mainly probed by f gas .We thus argue that the gas fraction can be used to set constraints on the cosmological model, albeit provided that one properly models the gas effects taking place inside clusters and provided that f gas is used in combination with other probes.MCMC results for a 10% shift of the depletion factor in our analysis in the VB scenario.We note that in this case, we assume a prior on the total value of the bias.

Fig. 1 .
Fig. 1.Gas mass fraction of Planck-ESZ clusters with respect to redshift (top) and with respect to cluster mass (bottom).The yellow bands in both plots mark the Planck Collaboration et al. (2020) value of Ω b /Ω m = 0.156 ± 0.003.

Fig. 2 .
Fig. 2. Binned mass-redshift plane of the Planck-ESZ sample.Inside each bin we compute the mean value of f gas .We show the number of clusters included in each bin and mark the delimitation of each subsample.

Fig. 5 .
Fig. 5. 1D and 2D posterior distributions for the bias parameters in the three mass-and redshift-selected samples.The levels of the contours mark the 68% and 95% confidence levels.Gray dashed lines mark the reference values (α, β) = (0, 0).
Fig.6.Fits obtained from our analysis, for the different mass-and redshift-selected samples.The results in the top two panels are represented with respect to redshift at constant mass (respectively minimal and maximal masses of the full sample), while the bottom two panels are the results presented with respect to mass, at a fixed redshift (respectively minimal and maximal redshifts of the full sample).The shaded areas around the curves mark the 68% and 95% confidence levels.The blue dashed line marks the reference value Ω b /Ω m = 0.156.

Fig. 7 .
Fig. 7. Comparison of the total masses inside the Planck-ESZ sample between XMM-Newton and Chandra.

5. 1 .
Possible sources of systematic effects 5.1.1.Instrumental calibration effects All of the masses used in this work were taken from Lovisari et al. (2020b), where the authors used XMM-Newton observations.Possible calibration effects may have impacted the mass measurements in their work and induced biases (see e.g.Mahdavi et al. (2013) or Schellenberger et al. (2015)) This effect has been taken into account under the form of a K = 1.0 ± 0.1 parameter in our analysis, but we try here to check if this assumption is sound.The clusters of the Planck-ESZ sample have also been observed using Chandra, in Andrade-Santos et al.

Fig. 9 .
Fig. 9. 1D and 2D posteriors obtained when comparing the bias parameters derived for the three samples, when considering a linear evolution of the bias.The contours mark the 68% and 95% c.l.
Fig. 10.Comparison between the PSZ2Cosmo sample used in Salvati et al. (2019) and the Planck-ESZ sample.
) meaning Ω b /Ω m ∝ B. When assuming a low value of the bias B = 0.62 ± 0.03 from Planck Collaboration et al. (2020), we show that the derived baryon fraction becomes Ω b /Ω m = 0.108±0.018,which is incompatible with the value of the baryon fraction derived from the CMB, Ω b /Ω m = 0.156 ± 0.003 (Planck Collaboration et al. 2020), as we show in Figure12.Assuming a low value of the bias would then imply that the Universe hosts roughly 20% less baryons than expected.We thus argue that a value of the bias B = 0.62 ± 0.03 seems highly unlikely provided our gas fraction data.

Fig. 12 .
Fig. 12.Comparison of the expected baryon fraction Ω b /Ω m between the bias we derived in the VB + Ω m scenario and B = 0.62 from (Planck Collaboration et al. 2020).
Fig. C.1.MCMC results for a 10% shift of the depletion factor in our analysis in the VB scenario.We note that in this case, we assume a prior on the total value of the bias.

Table 1 .
Set of priors used in our analysis.

Table 4 .
Constraints obtained on bias and cosmological parameters, when assuming a linear parametrization.Uncertainties are given at 68%c.l.