1 Introduction

Galaxy clusters are the most massive virialized objects in the universe and provide wonderful laboratories for studying a wide range of topics from galaxy evolution to cosmology [1, 2] to fundamental Physics [3]. They are also key in pinning down the dark energy equation of state, as they probe both the geometry of the Universe and the growth of structure [4, 5]. Galaxy clusters are therefore, one of the flagship probes for Stage IV dark energy experiments [6]. In the past two decades, a large number of galaxy clusters have been discovered upto very high redshifts, thanks to myriad multi-wavelength surveys, which have greatly facilitated these science goals.

Since galaxy clusters are the largest gravitationally bound objects, their matter content present in the intracluster medium has provided a powerful probe to constrain the cosmological parameters. The X-ray emitting gas in the intra-cluster medium dominates the baryon budget in clusters. The idea of using the gas mass fraction as a cosmological tool was first introduced by White et al. [7]. The gas mass fraction, which is defined as \(f_{gas} \equiv M_{gas} / M_{tot}\) is expected to match closely the cosmic baryon fraction, viz. \(\Omega _b / \Omega _m\), where \(\Omega _b\) and \(\Omega _m\) are the cosmic baryon density and total matter density respectively, in units of critical density (\(\rho _c\)). There have been a large number of studies, constraining the matter density (\(\Omega _m\)), dark energy equation of state parameter (\(\omega \)), and the curvature density (\(\Omega _k\)) using the \(f_{gas}\) measurements [8,9,10,11,12,13,14,15,16]. The \(f_{gas}\) test has also been used to test alternatives to \(\Lambda \)CDM such as \(R_h=ct\) universe [17]. The main assumption in these studies is that the gas fraction doesn’t evolve with redshift for hot, massive and relatively relaxed clusters. However, the evolution of \(f_{gas}\) as a function of redshift needs to be determined in a model-independent method, in order to probe the robustness of these results. This evolution is usually parametrized using the gas depletion factor \(\gamma \), which is the ratio by which the baryonic gas fraction in galaxy cluster gets depleted by the cosmic mean baryon fraction, and is given by: \(\gamma = f_{gas} (\Omega _b / \Omega _m)^{-1}\) [13]. So the study of gas depletion factor and it’s evolution with redshift play a pivotal role in the robustness of the \(f_{gas}\) test. As argued in [18], the results for the gas depletion factor obtained using only the data could then be used as a prior for any cosmology studies with \(f_{gas}\).

Therefore, a large number of studies have been undertaken using both simulations and data to test the variation of the gas depletion factor \(\gamma (z)\) as a function of redshift. This evolution is usually modeled using the following function.

$$\begin{aligned} \gamma (z)= \gamma _0(1+\gamma _1 z) \end{aligned}$$
(1)

where \(\gamma _0\) is the normalization factor and \(\gamma _1\) represents the evolution of gas depletion factor with redshift.

Battaglia et al. [19] studied the measurement biases of \(f_{gas}\) using hydrodynamical simulations including different kinds of processes like star formation, AGN feedback and radiative physics. They found a constant gas mass fraction at \(R_{500}.\)Footnote 1 Thereafter, Planelles et al. [21] studied the baryon fraction within \(R_{2500}\), \(R_{500}\), and \(R_{200}\) by using three different suites of simulations of galaxy clusters in the redshift range \(0 \leqslant z \leqslant 1\) for different physical processes in clusters. The first set of simulations (NR) involved non-radiative hydrodynamical simulations. The second suite (CSF) studied the effects of cooling, supernovae feedback, and star formation. Finally the third set of simulations (AGN) included the same physics as CSF, along with AGN feedback. Furthermore, they discussed the baryon and gas depletion factor and its dependence on the cluster radius, redshift and baryonic physics. They concluded that the evolution of the depletion factor is negligible with redshift (for \(z < 1\)) regardless of the cluster radius and baryonic physics. This aforementioned work obtained: \(\gamma _0 = 0.79\pm 0.07\) and \(\gamma _1 = 0.07\pm 0.12\) at \(R_{2500}\), \(\gamma _0 = 0.85\pm 0.03\) and \(\gamma _1 = 0.02\pm 0.05\) at \(R_{500}\), \(\gamma _0 = 0.86\pm 0.02\) and \(\gamma _1 = 0.00\pm 0.04\) at \(R_{200}\) [21]. We note that all these aforementioned simulations assumed \(\Lambda \)CDM as the base cosmological model.

Holanda et al. [18] made the first attempt to study the gas depletion factor using only observations. They carried out an analysis using 40 \(f_{gas}\) measurements [16] observed by the Chandra X-ray telescope spanning the redshift range \(0.078 \leqslant z \leqslant 1.063\) characterized as massive, morphologically relaxed systems and with \(kT \geqslant 5\) keV. In their work, \(f_{gas}\) measurements were obtained from inside a shell with radii between (0.8–1.2) \(R_{2500}\). Their analysis assumed the validity of cosmic distance duality relation (CDDR) [22], and the use of different distance indicators, e.g. Type Ia SNe (from Union 2.1 and JLA compilation) and \(\Lambda \)CDM model. They obtained \(\gamma _0 = 0.86\pm 0.04\) and \(\gamma _1 = -0.04\pm 0.12\) [18]. Therefore, their analysis showed no evolution of the gas depletion factor with redshift.

Then, in a followup work, Holanda [23] used 38 Chandra X-ray \(f_{gas}\) measurements from [24] in the redshift range \(0.14 \leqslant z \leqslant 0.89\), along with angular diameter distances from X-ray/SZ measurements. Unlike [18], they did not use CDDR to derive the angular diameter distance. They reported \(\gamma _0\) = \(0.76\pm 0.14\) and \(\gamma _1\) = \(-0.42_{-0.40}^{+0.42}\) at \(R_{2500}\), which also indicates no time evolution for the gas depletion factor.

Most recently, Zheng et al. [25] (Z19, hereafter) did a similar study using data for 182 clusters with the SZ [26,27,28] selected sample from the Atacama Cosmology Telescope Polarization experiment (ACTPol) spanning the redshift range \(0.1 \leqslant z \leqslant 1.4\) [29]. The main difference between Z19 and the previous works [18, 23] is that Z19 considered the data at \(R_{500}\), whereas the latter considered \(R_{2500}\). However, Z19 did not use direct \(f_{gas}\) measurements, but instead resorted to a semi-empirical relation between the gas mass fraction and the total cluster mass [30]. To estimate the angular diameter distance (needed for the estimation of the gas depletion factor), a model-independent approach using Gaussian Processes regression was used. They reported that the gas depletion factor decreases as a function of redshift. Similar to Z19, we study the evolution of the gas depletion factor for two different SZ selected cluster samples. The first set comprises a sample 94 South Pole Telescope (SPT) selected clusters in the redshift range \(0.278 \leqslant z \leqslant 1.32\), and the second set uses 120 Planck Early SZ (ESZ) data covering the redshift range \(0.059 \leqslant z \leqslant 0.546\).

Fig. 1
figure 1

\(f_{gas}\) data of 94 SPT-SZ clusters sample taken from [31, 32] and 120 Planck ESZ galaxy clusters sample taken from [46] as a function of z

Full size image

This paper is structured as follows. In Sect. 2, we briefly describe both our cluster samples, and their \(f_{gas}\) measurements. In Sect. 3, we describe the procedure used to determine the evolution of gas mass fraction. Our analysis along with results are discussed in Sect. 4. A discussion of our results can be found in Sect. 5. Finally, we conclude in Sect. 6.

2 Galaxy cluster sample and their \({{\varvec{f}}_{gas}}\) measurements

2.1 For SPT-SZ sample

The SPT-SZ cluster sample used for this analysis, consists of 91 SPT clusters from [31] (C18, hereafter) and three additional SPT-SZ clusters: namely SPT-CLJ 0205-5829, SPT-CLJ 0615-5746 and SPT-CLJ2040-5726 from [32], giving a total of 94 clusters, with a mass threshold of \(M_{500} \geqslant 2.5 \times 10^{14}\) \(M_\odot \), and in the redshift range of \(0.278 \leqslant z \leqslant 1.32\). These clusters were detected in the 2500 deg\(^2\) South Pole Telescope (SPT) SZ survey [33]. In this work, we used redshifts for each of these clusters from Bocquet et al. [34]. The SPT is a 10 m telescope at the South Pole, that has imaged the sky at three different frequencies, viz. 95 GHz, 150 GHz, and 220 GHz [33]. The SPT collaboration carried out a 2500 square degree survey between 2007 and 2011 to detect galaxy clusters using the SZ effect. This SPT-SZ survey detected 516 galaxy clusters with a mass threshold of \(3 \times 10^{14} M_{\odot }\) upto redshift of 1.8 [34, 35]. Detailed properties of the SPT clusters are discussed in [35]. Their redshifts have been obtained using survey data as well as pointed spectroscopic and photometric observations [36,37,38,39,40]. The first step in estimating the gas mass fraction consists of determining the total mass of the cluster at \(R_{500}\), known as \(M_{500}\). This was estimated from the SZ detection significance or signal to noise ratio, using SZ observable to mass scaling relations discussed in [41, 42]. These scaling relations have been self-consistently calibrated using a combination of X-ray, weak lensing, and number count distributions [31]. Here, we briefly describe the procedure used to derive the \(f_{gas}\) measurements below. More details can be found in the aforementioned works.

C18 estimated the total mass (denoted by \(M_{500}\)) and gas mass (which they refer to as intra-cluster medium mass, denoted by \(M_{ICM}\)) from the analysis of Chandra X-ray observations. All the 91 clusters have been imaged with Chandra, either through the Chandra X-ray Visionary Project led by the SPT collaboration as well as other proposals led by non-SPT members.

To estimate \(M_{ICM}\), a modified-\(\beta \) [43] model is used to fit the gas density, which was estimated from the observed surface brightness profile in the energy range 0.7–2.0 keV out to 1.5 \(R_{500}\). This surface profile was then fit to the modified \(\beta \)-model [43]. More details of this analysis can be found in [44]. We then calculate the gas mass fraction \(f_{gas}\), as the ratio of \(M_{ICM}\) to \(M_{500}\). This estimated gas mass fraction as a function of redshift is shown in Fig. 1. We note that C18 also investigated the trends of gas mass fraction as a function of redshift. However the parametric form which they have considered (Eq. 6 in C18) is different than that considered in this work (cf. Eq. (5)), and their analysis also includes a dependence on \(M_{500}\). A comparison between our results will be discussed in Sect. 5.

Another thing to point out is that for calculating this gas mass fraction, \(M_{500}\) was estimated using the observable SZ to mass scaling relation derived in a previous SPT work [42]. However, this scaling relation can change when using the cosmology priors from the Planck CMB anisotropy observations [41]. If this cosmology is adopted, \(M_{500}\) would be larger by about 22%, which in turn would also effect the baryon fraction numbers computed in C18. However, the Planck Cosmology prior is disfavored, since the weak lensing and dynamical masses are in agreement with the standard SPT cluster analysis [31]. In this work, we therefore do not consider the impact from the Planck cosmology prior.

Table 1 Redshifts of Planck-ESZ clusters used for this analysis
Full size table

2.2 For Planck ESZ sample

In this work, we use the gas mass fraction measurements of 120 Planck Early SZ (ESZ) clusters [45] spanning the redshift range \(0.059 \leqslant z \leqslant 0.546\) [46]. All these clusters were also imaged with XMM-Newton up to \(R_{500}\) [47]. The only difference compared to this dataset is that we used updated redshifts from Planck SZ2 catalog [48], instead of the ESZ redshifts used in [47]. These redshifts have also been obtained using pointed optical follow-ups with various telescopes (e.g. [49]) and are still been continuously refined (e.g. [50, 51]). Therefore, the redshifts of some of these clusters have further been updated after the publication of the Planck SZ2 catalog. The redshifts used for this analysis can be found in Table 1.

Lovisari et al. [46] derived the total mass of the cluster at \(R_{500}\) from the equation of hydrostatic equilibrium and assuming spherical symmetry:

$$\begin{aligned} M(r<r_{500}) = -\frac{kT_{gas}}{G\mu m_p{}{}} {\left[ \frac{d\log \rho _{gas}}{d\log r} + \frac{d\log T_{gas}}{d\log r}\right] } \end{aligned}$$
(2)

where \(T_{gas}(r)\) and \(\rho _{gas}(r)\) denote the temperature and density profiles (see [46, 47]), k is the Boltzmann’s constant and \(\mu m_p\) represents the mean molecular weight of the intra-cluster gas.

The gas density was obtained from the X-ray surface brightness profile in the 0.3–2 keV energy band after fitting a double-\(\beta \) model [47]. The temperatures were obtained by fitting the X-ray spectra to a thermal plasma emission model [47]. The gas mass \(M_{gas}\) was then estimated by integrating the density profile upto \(R_{500}\) after assuming spherical symmetry. Both \(M_{gas}\) and \(M_{500}\), along with their errors are tabulated in Lovisari et al. [46]. Thereafter, their ratio \(M_{gas}/M_{500}\) gives the desired gas mass fraction \(f_{gas}\) values for each cluster. Figure 1 shows the derived \(f_{gas}\) values for Planck ESZ clusters as a function of redshift.

3 Method

The gas mass fraction in galaxy clusters has been used as a cosmological probe. The expression for \(f_{gas}\) with respect to a reference cosmological model is given by [13]:

$$\begin{aligned} f_{gas}^{ref}= K(z)A(z)\gamma (z)\left( \frac{\Omega _b}{\Omega _m}\right) {\left[ \frac{D_A^{ref}(z)}{D_A(z)}\right] }^{3/2} \end{aligned}$$
(3)

where \(\Omega _b\) and \(\Omega _m\) denote the cosmic baryon and matter fraction density (obtained from Planck 18 Cosmological analysis [58]) scaled by the critical density \(\rho _c\); K(z) = \(1.0\pm 0.1\) [13] is the instrument calibration constant, which also accounts for any bias in mass due to non-thermal pressure and bulk motions in baryonic gas; A(z) represents the angular correction factor which is almost in all cases, close to unity, and hence can be neglected. The terms in the square parenthesis denote the variation in \(f_{gas}\) value as the underlying cosmology is changed. \(D_{A}\) is the angular diameter distance and ref indicates the fiducial cosmology. For our reference cosmology, we have used a flat \(\Lambda \)CDM cosmology (\(\Omega _m\) = 0.3 and \(H_0 = 70\) km/s/Mpc [46]).

Assuming flat \(\Lambda \)CDM cosmology, the angular diameter distance for a reference cosmology can be calculated as follows [59]:

$$\begin{aligned} D_{A}^{ref} = \frac{c}{H_0} \left( \frac{1}{1+z}\right) \int _{0}^{z}\frac{dz^{'}}{\sqrt{\Omega _m{(1+z^{'})^3}+(1 - \Omega _m)}} \end{aligned}$$
(4)

From Eq. (3), we can write the gas depletion factor, \(\gamma (z)\) as follows:

$$\begin{aligned} \gamma (z)= \frac{f_{gas}^{ref}}{ A(z) K(z)} \left( \frac{\Omega _m}{\Omega _b}\right) {\left[ \frac{D_A(z)}{D_A^{ref}(z)}\right] }^{3/2} \end{aligned}$$
(5)

This equation shows that \(\gamma (z)\) is sensitive to measurements of angular diameter distance. In this work, similar to Z19 we have used Gaussian Process Regression to derive the angular diameter distances by using H(z) data from cosmic chronometers, which are agnostic to the underlying cosmology. \(D_A(z)\) is given by

$$\begin{aligned} D_{A} (z) = \left( \frac{c}{1+z}\right) \int _{0}^{z}\frac{dz^{'}}{H(z')} \end{aligned}$$
(6)

where \(H(z')\) is the non-parametric estimate of the Hubble parameter, whose determination will be described in the next section.

Fig. 2
figure 2

Reconstruction of H(z) using Gaussian Processes Regression. The data points represent the 31 H(z) cosmic chronometer measurements taken from Li et al. [60]. The blue line indicates the best GP fit to data along with \(1\sigma \) and \(2\sigma \) error bands shown by two different shades of green bands

Full size image

4 Analysis and results

Cosmic chronometers provide a model-independent measurement of H(z), based on the age difference between two passively evolving galaxies [61]. Therefore, they have been widely used for a variety of tests of the standard cosmological model, as well as a whole suite of cosmological measurements [62,63,64,65,66,67,68,69,70,71,72]. For this work, we used cosmic chronometer data to provide a non-parametric estimate of the expansion history, needed to evaluate Eq. (6).

Gaussian Processes (GPs) extend the idea of Gaussian distribution, and are characterized by the mean function and the covariance (kernel) function. They offer a non-parametric way to model a function [62]. In this work, we choose a squared exponential (RBF) kernel function which is defined as,

$$\begin{aligned} K(x,\tilde{x}) = \sigma _f^2 exp {\left[ \frac{-(x-\tilde{x})^2}{2l^2}\right] }, \end{aligned}$$
(7)

where \(\sigma _f^2\) and l are the hyperparameters of the kernel function. The length parameter l controls the smoothness of the function. To implement the GPs, we used the GaPP (Gaussian Process in Python) code [73].

Following Z19, we used the 31 H(z) cosmic chronometer measurements from  Li et al. [60] (same as that used in [72]) in the redshift range \(0.07 \leqslant z \leqslant 1.965\) to obtain non-parametric estimates of H(z) at any redshift. Figure 2 shows the GP reconstructed H(z) from cosmic chronometers along with \(1\sigma \) and \(2\sigma \) uncertainties. We then reconstruct the GP reconstructed H(z) values for SPT-SZ and Planck ESZ redshifts in order to obtain the angular diameter distances. Thereafter, these H(z) values were used to derive the angular diameter distances \(D_A\) using Eq. (6). \(D_A^{ref}\) is calculated from Eq. (4) by assuming the cosmological parameters outlined in Sect. 3.

Finally, we reconstructed the gas depletion factor \(\gamma (z)\) using Gaussian Process, after incorporating all the other factors presented in Eq. (5). We also propagated the errors in \(\Omega _b\), \(\Omega _M\), K, angular diameter distance, and redshifts, if they were provided. Figures 3 and  4 display our reconstructed \(\gamma (z)\) as a function of z for both SPT-SZ and Planck ESZ data, respectively. Both these data show opposite trends for the gas depletion factor with redshift. We find a decreasing \(\gamma \) for SPT-SZ data, whereas for Planck ESZ data, a slightly increasing trend is observed. For comparison, Planelles et al. [21] reported a negligible evolution of gas depletion factor by using Eq. (1) for their NR simulations. They obtained \(\gamma _0 = 0.85\pm 0.03\) and \(\gamma _1 = 0.02\pm 0.05\) at \(R_{500}\). The simulation results can be found in the red shaded region in Figs. 3 and  4. Therefore, at face value, our results do not agree with these hydrodynamical simulations based on \(\Lambda \)CDM.

Fig. 3
figure 3

Reconstructed \(\gamma (z)\) as a function of z. The blue line represents the best GP reconstructed fit for the SPT-SZ data along with the \(1\sigma \) error band (shown by blue shaded region). The red band shows the NR simulation results at \(R_{500}\) from [21] by assuming the Eq. (1)

Full size image
Fig. 4
figure 4

Reconstructed \(\gamma (z)\) as a function of z. The green line represents the best GP reconstructed fit for the Planck ESZ data along with the \(1\sigma \) error band(shown by green shaded region). The red region shows the NR simulation results at \(R_{500}\) from [21] by assuming the Eq. (1)

Full size image

To confirm and quantify the significance of the trends seen in \(\gamma (z)\), we also did a parametric fit of our estimated depletion factor measurements to Eq. (1). We maximize the following likelihood function \({\mathcal {L}}\) in our analysis, where

$$\begin{aligned} -2\ln {\mathcal {L}} = \sum _{obs} \ln 2\pi {\sigma _{obs}^2}+ \sum _{i=1}^{n}\frac{(\gamma _{model}(z_i)-\gamma _{obs}(z_i))^2}{\sigma _{i,obs}^2} ,\nonumber \\ \end{aligned}$$
(8)

where \(\gamma _{model}\) is defined in Eq. (1), \(\gamma _{obs}\) is obtained from Eq. (5), and \(\sigma _{i,obs}^2\) represents the error in \(\gamma _{obs}\). \(\sigma _{i,obs}\) is calculated by propagating the errors in \(f_{gas}^{ref}\), K, \(\Omega _m\), cluster z, \(\Omega _b\), \(D_A(z)\), and \(D_A^{ref}\). In addition to these observational errors, \(\sigma _{obs}\) also includes an unknown intrinsic scatter (\(\sigma _{int}\)) added in quadrature, similar to our recent works on galaxy cluster scaling relations [74,75,76].

To estimate the model parameters(\(\gamma _0\) and \(\gamma _1\)), we maximize the likelihood using the emcee MCMC sampler [77]. For both the datasets, we adopt uniform priors on \(\gamma _0\) and \(\gamma _1\): \(-1.5 \leqslant \gamma _0 \leqslant 1.5\) and \(-2.0 \leqslant \gamma _1 \leqslant 3.0\). For the intrinsic scatter, similar to [78], we assume log-uniform priors between \(10^{-5}\) and 0.1. The log-uniform prior ensures that the intrinsic scatter is greater than zero. The 68%, 95%, and 99% 2-D marginalized credible intervals, along with the marginalized one-dimensional likelihoods for each parameters are displayed in Figs. 5 and 6 respectively.

Fig. 5
figure 5

For SPT-SZ dataset: constraints on the parameters \(\gamma _0\) and \(\gamma _1\) along with \(\ln (\sigma _{int})\), defined in Eq. (8). The plots along the diagonal are the one-dimensional marginalized likelihood distributions. The contour plots represents the two-dimensional marginalized constraints showing the 68%, 95%, and 99% credible regions. These contours have been obtained using the Corner python module [79]

Full size image
Fig. 6
figure 6

For Planck ESZ dataset: constraints on the parameters \(\gamma _0\) and \(\gamma _1\) along with \(\ln (\sigma _{int})\). The plots along the diagonal are the one-dimensional marginalized likelihood distributions. The contour plot represents the two-dimensional marginalized constraints showing the 68%, 95%, and 99% credible regions. These contours have been obtained using the Corner python module [79]

Full size image

From this analysis, we get \(\gamma _0 = 0.915\pm 0.050\) and \(\gamma _1 = -0.305_{-0.058}^{+0.064}\) for SPT-SZ data, and \(\gamma _0 = 0.741\pm 0.029\) and \(\gamma _1 = 0.834_{-0.199}^{+0.214}\) for Planck ESZ data. The SPT-SZ data results indicate a decreasing trend of \(\gamma \)(z) with redshift at about 5\(\sigma \) significance. On the other hand, the Planck ESZ results show an increasing \(\gamma \) as a function of redshift at about \(4\sigma \). Both the results contradict those from hydrodynamical simulations [21]. Table 2 summarizes the constraints on parameters \(\gamma _0\) and \(\gamma _1\), which we have obtained from the datasets, along with a summary of previous compilations in literature using both data and simulations. We note that among the previous observational results, only Z19 had found \(\gamma (z)\) decreasing with redshift at about \(1.6\sigma \).

Table 2 Constraints on parameters \(\gamma _0\) and \(\gamma _1\) (cf. Eq. (1)) from different estimates in literature (including this work). The last two rows show the gas depletion factors for a subset of SPT-SZ and Planck ESZ clusters within the same redshift range
Full size table

5 Discussions

Given the contradictory results for \(\gamma (z)\) between the SPT-SZ and Planck-ESZ datasets as well as with previous results in literature, we carry out various sanity checks of our results along with extensions of tests carried out earlier, to verify the robustness of our conclusions, and to see if there is one deciding parameter, which governs the difference between the results.

5.1 Comparison with C18

A test of the variation of the gas mass fraction with redshift (using different parametric forms) for the SPT-SZ sample was also carried out in C18, and no evidence for a variation of the depletion factor for redshift was found (cf. Table 6 in C18.) At prima facie, it may seem that our results contradict those in C18. However, \(\gamma (z)\) defined in this work (cf. Eq. (1)) is a function of \(\Omega _b\), \(\Omega _M\), and \(D_A(z)\), besides the gas mass fraction. On the other hand, C18 carried out a joint test for the variation of gas mass fraction with redshift and \(M_{500}\). The functional forms used in C18 to test for redshift trends are \(\propto (1+z)^{\alpha }\) and \((E(z))^{\gamma }\), where E(z) is the cosmic expansion factor in \(\Lambda \)CDM [70]. Here, \(\gamma (z)\) estimated in this work is proportional to \(\left[ D_A(z)\right] ^{3/2}\) (cf. Eq. (5)), which is the integral of the reciprocal of E(z). \(\gamma (z)\) studies in this work also does not have an explicit dependence on the mass. We have also reconstructed \(D_A(z)\) in a non-parameteric way using cosmic chronometers.

However, in order to test the consistency of our results with C18 (for SPT-SZ), we now check for a variation of only the raw \(f_{gas}\) data to a few different parametric functions (including those used in C18). The full list of functions used to test variations of gas mass fraction are indicated below:

$$\begin{aligned} f_{gas}= & {} f_0(1 + f_1 z) \end{aligned}$$
(9)
$$\begin{aligned} f_{gas}= & {} \left( \frac{1+z}{1+z_{pivot}}\right) ^{\alpha } \end{aligned}$$
(10)
$$\begin{aligned} f_{gas}= & {} \left( \frac{E(z)}{E(z_{pivot})}\right) ^{\beta } \end{aligned}$$
(11)

where \(z_{pivot}\) is the pivot value of the redshift, equal to 0.6 for SPT and 0.2 for Planck. With these functional forms, we remove any explicit dependence on the angular diameter distance ratio. The best-fit values of the free parameters in the above equations for both the datasets were found by constructing a likelihood similar to Eq. (8). These best-fit parameters can be found in Table 3 for both SPT-SZ and Planck-ESZ. We find that for the SPT-SZ sample, the parameters \(f_1\), \(\alpha \), and \(\beta \), which encode the dependence on redshift are consistent with no variation (within 1\(\sigma \)). Therefore, this agrees with the results in C18, which found no redshift dependence.

However, for the Planck sample \(f_1\) in Eq. (9) is equal to \(0.85 \pm 0.17\), indicating a positive slope which is 5\(\sigma \) discrepant from zero (which corresponds to no redshift evolution). This is one of the main determining factors for the positive \(\gamma _1\) for Planck-ESZ. The deviations for \(\alpha \) and \(\beta \) from a zero value are within \(1-2\sigma \).

5.2 \(\gamma (z)\) within the overlapping redshift range

The SPT sample extends upto very high redshift (\(z=1.32\)) compared to Planck, whose most distant cluster is located at \(z=0.546\). Similarly, the Planck sample contains about 84 low redshift clusters with \(z<0.278\), which is the lowest redshift for the SPT sample. We now check if the gas depletion factor for both the samples within the overlapping redshift range (\(0.278\leqslant z \leqslant 0.546\)) display the same trend. We get a total of 34 SPT and 36 Planck clusters with these redshift cuts. We applied the same procedure as in Sect. 4 to determine \(\gamma _0\) and \(\gamma _1\) for these sub-samples. These values can be found in the last two rows of Table 2. The best-fit value of \(\gamma _1\) for SPT-SZ and Planck-SZ is equal to \(1.02_{-0.79}^{+1.44}\) and \(0.24_{-0.47}^{+0.7}\), respectively. Therefore, we see that although there is a mild deviation of \(\gamma _1\) from a zero value for SPT-SZ, their values are consistent between the two datasets within \(1\sigma \).

5.3 Dependence on \(M_{500}\)

The fits to \(f_{gas}\) done in C18 also had an explicit dependence on \(M_{500}\) given by \(f_{gas} \propto M_{500}^{\zeta }\). C18 then found \(\zeta =0.33 \pm 0.07\) with a \(4.5\sigma \) deviation from self-similarity (where, self-similarity corresponds to \(\zeta =0\)) in agreement with previous works [30, 43, 80, 81]. So far, we have not accounted for any dependence on the halo mass, while estimating the gas depletion factor as a function of redshift.

Table 3 Constraints on the parameters of various extended tests of the gas depletion factors as discussed in Sect. 5 for both SPT and Planck clusters
Full size table

In order to check the possible dependence on the mass of the cluster, and evaluate its impact on the gas depletion factor, we augment Eq. (1) with an additional dependence on \(M_{500}\) as follows:

$$\begin{aligned} \gamma (z)= \gamma _0(1+\gamma _1 z)\left( \frac{M_{500}}{M_{pivot}}\right) ^{\zeta } \end{aligned}$$
(12)

where \(M_{pivot}\) corresponds to the median mass for the SPT and Planck samples. These are equal to \( 4.8 \times 10^{14}M_{\odot }\) and \(6.0 \times 10^{14}M_{\odot }\) for SPT-SZ and Planck ESZ, respectively. To determine these parameters, we use an extension of the same procedure as in Sect. 4, by including an additional dependence on \(M_{500}\). While constructing the likelihood, the errors in \(M_{500}\) were also included. The marginalized posterior intervals for the gas depletion parameters along with \(\zeta \), which encodes the variation with \(M_{500}\), can be found in Figs. 7 and 8 for SPT-SZ and Planck-ESZ, respectively. The best-fit values for these parameters can be found in Table 3. For SPT-SZ, we find that \(\zeta = 0.257 \pm 0.09\), indicating a 2.9\(\sigma \) deviation from self-similarity. For Planck-ESZ, we have \(\zeta =0.083 \pm 0.06\), with only a 1.4\(\sigma \), deviation from self-similarity. After allowing for this dependence on \(M_{500}\), we find that \(\gamma _1=-0.23 \pm 0.08\) and \(0.57 \pm 0.26\) for SPT-SZ and Planck-ESZ respectively. This shows that the discrepancy with respect to a constant depletion factor as a function of redshift decreases to \(2.9\sigma \) and \(2.2\sigma \) for SPT-SZ and Planck-ESZ, respectively, as compared to the \(4-5\sigma \) significance, which we had obtained earlier without including an explicit dependence on the halo mass. Therefore, we conclude that although the statistical significance of the rising/falling slope as a function of redshift decreases, it does not completely go away, when we include a dependence on the total mass. Only its significance reduces to between \(2-3\sigma \). However, given that the inclusion of halo mass in our fit reduces the discrepancy with respect to no evolution of \(\gamma _1\), one cannot decouple the two factors while testing for a redshift of the gas depletion factor.

5.4 Other possible sources of systematics

Here, we list other possible sources of errors, which could affect our estimates of \(\gamma (z)\). A detailed study of the impact of each of these effects is beyond the scope of this work. The most crucial ingredient is the determination of cluster masses, which are obtained from the observable to halo mass scaling relations and is limited by the systematic uncertainties. The gas mass can be sensitive to the assumed temperature profile [15]. Errors due to hydrostatic equilibrium assumption could be upto 15–20% [82]. Other possible sources of systematic errors in the gas mass fraction determination such as magnetic field, thermal evaporation, and non-linear translation of analysis variables are discussed in [83,84,85]. Finally, we point out that the list of SPT-SZ and Planck-SZ analyzed in this work is not the full sample. The PSZ2 catalog contains 1653 cluster candidates (with 1,203 confirmed clusters) with signal to noise ratio above 4.5 [48]. This work has analyzed only about 10% of the full sample, which was followed up in X-rays at the time of Planck ESZ release. Optical and X-ray studies of Planck SZ clusters is still in progress [51, 86]. Similarly, the SPT-SZ dataset analyzed in this work is only a quarter of the full 2500 sq. degree survey sample used for the cosmological analysis in [34]. However, the SPT clusters are also been imaged in X-rays through dedicated Chandra [44] and XMM [87] based follow-ups. Since our currently analyzed subset may not be representative of the full SPT and Planck SZ sample, it would be interesting to carryout a followup analysis, once the follow-up data for the full sample is available.

6 Conclusions

The gas mass fraction (\(f_{gas}\)) in galaxy clusters has been used as a cosmological probe in a number of works. A key assumption in these analyses is that the gas mass fraction (at the particular radius used for cosmological measurements) does not vary with redshift. This has also been confirmed with hydrodynamical simulations of galaxy clusters which use \(\Lambda \)CDM as the base model, which show no evolution of the gas depletion factor with redshift (\(z< 1\)) [21].

Fig. 7
figure 7

Marginalized credible intervals (68%, 95%, 99%) for \(\gamma _0\), \(\gamma _1\), \(\zeta \), and \(\ln (\sigma _{int})\) for the SPT-SZ dataset. The parameter \(\zeta \) encodes the variation of the gas depletion factor as a function of \(M_{500}\) (cf. Eq. (12))

Full size image
Fig. 8
figure 8

Same as Fig. 7 for the Planck ESZ dataset

Full size image

To verify this ansatz with real data in a model-independent method without any dependence on the underlying cosmology, Z19 carried out a model-independent test for the variation of gas depletion factor using 182 clusters from ACTPol. Using the parameterization in Eq. (1), they found that \(\gamma _0 = 0.840\pm 0.025\) and \(\gamma _1 = -0.072_{-0.049}^{+0.044}\), which indicates a decreasing trend at about \(1.4\sigma \). In this work, we apply the same procedure as Z19 using the data sets from SPT-SZ and Planck ESZ in the redshift range \(0.278 \leqslant z \leqslant 1.32\) and \(0.059 \leqslant z \leqslant 0.546\), respectively. One difference with respect to Z19, however is that we used direct gas mass fraction measurements from SPT and Planck, instead of the empirical relation between halo and gas mass used in Z19. These measurements were obtained using a combination of SZ data along with follow-up X-ray observations from Chandra and XMM, for SPT and Planck, respectively.

Our \(f_{gas}\) measurements as a function of redshift can be found in Fig. 1 for SPT-SZ and Planck ESZ clusters. We derived the angular diameter distance for both the datasets using GP fitting from H(z) data reconstructed using cosmic chronometers. The reconstructed gas depletion factor can be found in Figs. 3 and 4 for SPT-SZ and Planck ESZ data respectively. To quantify the significance of our results obtained from GP reconstruction, we also did a maximum likelihood fit using Eq. (1). We find that \(\gamma _0 = 0.915\pm 0.05\) and \(\gamma _1 = -0.305_{-0.058}^{+0.064}\) for SPT-SZ data; whereas \(\gamma _0 = 0.741\pm 0.029\) and \(\gamma _1 = 0.834 \pm 0.2\) for Planck ESZ data. The credible intervals for \(\gamma _0\) as well as \(\gamma _1\) along with marginalized 1-D distributions are shown in Figs. 5 and 6 for SPT-SZ and Planck ESZ data respectively. We find \(\gamma \) decreases with redshift (at about 5\(\sigma \)) for SPT-SZ clusters, whereas it increases with redshift for Planck ESZ clusters (at about 4\(\sigma \)). Both these findings also contradict hydrodynamical simulations which show a constant gas depletion factor. Previously, only Z19 had found a slightly decreasing trend with redshift for \(\gamma (z)\) using ACT-Pol data.

Given the contradictory results between the two datasets, we carried out a series of cross-checks on our results using various extended tests, as outlined in Sect. 5. We used simple parametric forms to model only the gas mass fraction, without any explicit dependence on the angular diameter distance (cf. Eqs. (9), (10), (11).) We then evaluated the gas depletion factor for both the datasets within the overlapping redshift range. Finally, we also added an explicit dependence on \(M_{500}\) in the parameterization used for the gas depletion factor (cf. Eq. (12)). The results of these extended tests are tabulated in Table 3 and in Figs. 7 (SPT-SZ) and 8 (Planck ESZ). We find that the coefficients for the redshift dependent terms in Eqs. (9), (10), (11) are consistent with no variation within 1\(\sigma \). For Planck-ESZ, \(f_1\) in Eq. (9) is 5\(\sigma \) deviant with respect to zero, which is one of the main underlying causes for a positive value for \(\gamma _1\). If we analyze the evolution of the gas depletion factor for a subset of SPT-SZ and Planck ESZ clusters within the same redshift range, we find that \(\gamma _1\) is consistent for both the datasets within \(1\sigma \) (cf. last two rows in Table 2). Finally, when we allow for a variation of the gas depletion factor with halo mass, we find that the decreasing/increasing trends with redshift for SPT-SZ (Planck ESZ) still persist, albeit with reduced significance of about \(2-3\sigma \).

Therefore, our investigations on a model-independent test of the evolution of \(f_{gas}\) reveal that the redshift evolution of the gas depletion factor depends on the redshift range, and also changes depending on whether one includes/excludes a dependence on the halo mass. One possible reason for these changes is that the sample analyzed here is not the full SPT-SZ or the Planck SZ sample as the optical/X-ray followups is still ongoing. Another possibility is that the error in \(f_{gas}\) also depends on additional systematics such as observable to halo mass relation (which in turn is somewhat cosmology dependent), gas temperature and density profiles, hydrostatic equilibrium assumption, magnetic fields, thermal evaporation, which have not been completely accounted for, and could be redshift dependent. Nevertheless, given the observed deviation from self-similarity found for both the cluster samples (and in C18) and the significant change in the redshift dependence of \(f_{gas}\) without reference to any underlying cosmology, one must include the dependence in halo mass for modelling any evolution of \(f_{gas}\) with redshift.

We note that, our results do not affect any cosmological parameter estimations with \(f_{gas}\), since most of those results were obtained with \(f_{gas}\) estimated at \(R_{2500}\) [2, 11,12,13], whereas we used \(f_{gas}\) measurements at \(R_{500}\). Only a handful of works have used gas fraction measurements at \(R_{500}\) for cosmology studies [9, 15]. However, if the conflicting results for the two SZ datasets at \(R_{500}\) persist when extended to the full SZ dataset from both the telescopes, then it implies that one cannot use \(f_{gas}\) values at \(R_{500}\) as a stand-alone probe for any cosmological tests.