© C. H. A. Logan et al. 2022

1. Introduction

Galaxy clusters are the largest gravitationally bound objects in the Universe, and thus they are used to study the growth of cosmic structures and as cosmological probes. Tight constraints on cosmological parameters can be determined from measurements of the mass function of clusters, or their baryon fraction, as a function of redshift (see e.g. Allen et al. 2011, for a review). However, the accuracy of these constraints crucially depends on the accuracy with which cluster masses can be measured. For example, the exponential drop in the number density of clusters at high masses means that their mass function is a sensitive cosmological probe, but an error of ∼10% on the cluster mass can lead to a factor of two difference in the expected space density¹. The measurement of galaxy cluster masses is challenging, as most of the cluster mass is in the form of dark matter. Mass estimation techniques thus probe the masses of clusters indirectly via the effect of the gravitational potential of the cluster on its intracluster medium (ICM), its member galaxies, or images of background galaxies.

Optical data were used to give the first mass estimates of galaxy clusters. Zwicky (1933) measured the line-of-sight velocities of a number of cluster galaxies in the nearby Coma cluster, and, by computing the velocity dispersion, used the virial theorem to estimate a mass for the cluster. In the late 1960s and 1970s, the first X-ray satellites were launched, and it was in the 1990s that the first X-ray samples of hundreds of clusters were produced (Ebeling et al. 1998, 2001; Böhringer et al. 2000, 2001) from which X-ray hydrostatic masses could be calculated (e.g. Reiprich & Böhringer 2002). The launch of XMM-Newton and Chandra in 1999 has led to even more precise X-ray masses being measured for large samples of clusters (e.g. Martino et al. 2014). The first work using weak gravitational lensing (WL) to map the dark matter distribution, undertaken by Tyson et al. (1990), was followed by a number of other papers (e.g. Blandford et al. 1991; Miralda-Escude 1991; Kaiser 1992), and in the last 3 decades it has become increasingly common to use WL as a mass estimation technique (see Umetsu 2020 for a review).

For cosmological studies, the calibration of masses is critical to understand the biases that may be inherent to different mass estimation techniques. Cosmology requires large samples of clusters with masses which in general must be estimated from simple mass proxies that must be calibrated with more direct, higher fidelity mass measurements. Thus, the accuracy of the mass scale from those measurements is crucial, and biases in the mass scale can manifest as biases in the derived cosmological parameters. For example, there is currently tension between the cosmological constraints (Ω_m and σ₈) from the Planck cluster number counts (Planck Collaboration XXIV 2016) and the Planck cosmic microwave background (CMB) results (Planck Collaboration XIII 2016). The Planck cluster masses were calibrated using hydrostatic masses, and the tension in the results could be due to a negative bias (i.e. the true mass is underestimated) in the hydrostatic masses, often referred to as the hydrostatic bias.

More recent analyses have found lower tension between the values of σ₈ obtained from the Planck CMB and from cluster number counts measured with SZE or X-ray cluster surveys (e.g. Douspis et al. 2019; Bocquet et al. 2019; Planck Collaboration VI 2020; Garrel et al. 2022; Salvati et al. 2022). The improved agreement comes from updates to the CMB analysis and the mass calibration of the clusters, with all analyses favouring a non-negligible hydrostatic bias. However the magnitude of the hydrostatic bias is not agreed upon, and remains a critical point for cluster studies.

Some level of hydrostatic bias is expected a-priori, due to the assumption that the only pressure support in the ICM is the thermal pressure measured in X-rays. If any other sources of pressure (e.g. from turbulence, bulk motions or clumping of the ICM, or from cosmic rays) are also present, this leads to the X-ray hydrostatic mass method underestimating the true cluster mass (Nagai et al. 2007; Lau et al. 2009; Nelson et al. 2014).

The level of hydrostatic bias is expected to be higher in unrelaxed clusters compared with relaxed clusters, as unrelaxed systems are expected to have more non-thermal pressure support from bulk motions and turbulence of the cluster gas. However, all clusters, regardless of dynamical state, will experience some amount of non-thermal pressure support, as gas and substructures are always infalling onto the cluster as they undergo constant growth, leading to non-thermal pressure due to the residual motion of the ICM and the turbulence created. Constraining turbulence in the ICM requires high resolution X-ray spectrometry, and the only mission thus far capable of this has been Hitomi. Before the mission failed, Hitomi observed the Perseus cluster, finding that in the core of the Perseus cluster the pressure support from turbulence was 4% of the thermodynamic pressure Hitomi Collaboration (2016), which on its own would not lead to a large value of hydrostatic bias (though we note that the Perseus cluster is a dynamically relaxed cluster and, of course, not representative of the cluster population as a whole).

The hydrostatic bias may also be a function of radius, as the infalling gas may lead to stronger gas motions and therefore significant non-thermal pressure support at the outskirts of clusters. Bonamente et al. (2013) and Fusco-Femiano & Lapi (2018) both find evidence for significant non-thermal pressure at the cluster outskirts (for A1835 and A2142 respectively), and work using hydrodynamical simulations suggests that the non-thermal pressure due to the bulk gas motion increases with radius (Nagai et al. 2007; Lau et al. 2009; Vazza et al. 2009, 2018; Battaglia et al. 2012; Rasia et al. 2012). Specifically, Nagai et al. (2007) and Rasia et al. (2012) both find that at low (< 0.2 R₅₀₀) and high radius (> R₅₀₀) the bias increases². Nagai et al. (2007) find the bias to increase from < 10% at R₅₀₀ to ∼30% at twice this radius. Studies that have used observational data agree with this general trend (Siegel et al. 2018; Eckert et al. 2019).

To get a grasp on the biases of the cluster mass estimation methods, comparisons between methods themselves are undertaken, and simulated data are also used to compare the known mass from the simulation to the mass that is recovered if observational methods are applied to synthetic observations of the simulated cluster. Comparisons of X-ray hydrostatic masses with WL masses can yield insights into the magnitude of the hydrostatic bias, as WL masses are not sensitive to the state of the ICM (in contrast to X-ray hydrostatic masses). However, WL masses are sensitive to mass along the line of sight to the observed cluster, which can lead to scatter of ∼20−30% and positive bias of up to ∼20% (Hoekstra 2001, 2003; Hoekstra et al. 2011; Becker & Kravtsov 2011; Hwang et al. 2014), which needs to be understood and taken into account in these comparisons. A further complication is that the uncertainty on WL masses increases with decreasing cluster mass and decreasing cluster redshift (Hoekstra 2003), and even for a cluster with mass 5 × 10¹⁴ M_⊙ at a redshift of 0.3, the one sigma uncertainty can be as large as 40%. Some comparisons between WL masses and X-ray hydrostatic masses suggest that X-ray hydrostatic masses underestimate the true mass by ∼20−30% (von der Linden et al. 2014; Donahue et al. 2014; Sereno et al. 2015; Hoekstra et al. 2015). Work using hydrodynamical simulations also suggest a significant hydrostatic bias of ∼10−30% (Rasia et al. 2006, 2012; Nagai et al. 2007; Lau et al. 2009; Nelson et al. 2014). However Gruen et al. (2014), Israel et al. (2014), Applegate et al. (2016) and Smith et al. (2016) find no significant evidence for hydrostatic bias when comparing X-ray hydrostatic masses to WL masses.

Comparisons between X-ray and WL masses from observations also support the idea that disturbed clusters will have a larger hydrostatic bias than relaxed clusters (Zhang et al. 2010; Mahdavi et al. 2013). Biffi et al. (2016) also find a difference in hydrostatic bias between relaxed and disturbed clusters using simulations (though find no difference between cool core and non cool core systems).

Comparisons between X-ray and caustic masses are much rarer. Caustic masses (Diaferio & Geller 1997; Diaferio 1999) are a useful alternative to WL masses for testing the accuracy of hydrostatic masses. Like WL masses they are not sensitive to the state of the ICM, and biases in the caustic method are, in principle, well understood from simulations (Serra et al. 2011). One drawback of using caustic masses for comparison is that they have a large scatter of ∼30% (Serra et al. 2011), which comes predominantly from the assumption of spherical symmetry (and thus the viewing angle of the cluster). Thus, for a meaningful comparison, a large cluster sample is needed (several tens of clusters). Another drawback of the comparison is that the caustic mass method was conceived to measure the mass in the outer region of clusters, from 1 to 3 R₂₀₀, where the X-ray hydrostatic mass method cannot be applied because the X-ray surface brightness is too low to return a reliable signal at these radii. However, Serra et al. (2011) use simulations to estimate the bias in the caustic method at lower radii, such that a meaningful comparison between the two cluster mass estimators can still be made.

Diaferio et al. (2005) was the first work to compare X-ray and caustic masses, for a sample of three clusters, finding that for two clusters the mass profiles agreed well, and for the other cluster they did not, as the cluster was clearly out of equilibrium. Much larger cluster samples (from observational data) were used by Maughan et al. (2016), Andreon et al. (2017a) and Foëx et al. (2017) to infer caustic masses to be ∼20% lower (at R₅₀₀), ∼15% larger (at R₅₀₀) and ∼30% larger (at R₂₀₀) than the hydrostatic masses respectively. Armitage et al. (2019)³ use simulation data to compare caustic and X-ray hydrostatic masses (at R₂₀₀) and find a similar mass ratio to Foëx et al. (2017). Ettori et al. (2019), however, find their caustic masses to be underestimated by ∼40% (at R₂₀₀) compared to their hydrostatic masses (both computed using observational data). These contradictory results clearly show that the situation is complex and unsettled.

There have also been comparisons of caustic masses with weak lensing – Diaferio et al. (2005) do so for a sample of three clusters, and Geller et al. (2013) do so for a sample of 19 clusters. Both works find the mass methods to agree (independent of the dynamical state of the cluster) at around the virial radius to within ∼30%, consistent within the errors associated with each mass method.

In this paper, we aim to compare hydrostatic and caustic masses for the largest sample so far used for these purposes, in order to investigate the accuracy of the methods and provide new constraints on the hydrostatic bias that are independent of those based on WL masses. We improve on our earlier work (Maughan et al. 2016) by expanding the sample from 16 to 44 clusters. We also include a more comprehensive analysis of the systematics, and our main conclusions are based on a subset of 14 clusters for which the systematics are minimised.

In this paper we use the following notation to represent the masses of the clusters. M_X and M_C denote the masses measured with the hydrostatic and caustic methods, respectively. When these masses are specifically measured within R₅₀₀, we denote them as M_500, X and M_500, C. Unless otherwise stated, these are always measured within the R₅₀₀ determined from the hydrostatic analysis.

The structure of the paper is as follows. In Sect. 2 we discuss the sample selection and data preparation. Section 3 details the data analysis and modelling. In Sect. 4 we present initial results for the full sample, and then our main results from the subset found to have the most reliable mass estimates. We discuss our results in Sect. 5, including a review of the systematics affecting the mass estimates, and the implications for the hydrostatic bias, and summarise our findings in Sect. 6. Throughout this paper we assume a WMAP9 cosmology of H₀ = 69.7 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.718, and Ω_m = 0.282 (Hinshaw et al. 2013). All measurement uncertainties correspond to 68% confidence level, unless explicitly stated otherwise. When referring to base 10 logarithm we use ‘log’ or ‘log₁₀’, and use ‘ln’ if referring to the natural logarithm.

2. Cluster sample

For our analysis, we used the X-ray flux-limited subset of the Hectospec Cluster Survey (HeCS; Rines et al. 2013), a spectroscopic survey of X-ray selected clusters. The HeCS sample was constructed by matching clusters selected using the ROSAT All-Sky Survey (RASS; Voges et al. 1999) with the imaging footprint of the SDSS Data Release 6 (Adelman-McCarthy et al. 2008). The SDSS data were used to select candidate cluster member galaxies for spectroscopic follow-up, which was performed with MMT/Hectospec (Fabricant et al. 2005).

The X-ray flux limit applied to the HeCS sample in order to create the flux-limited HeCS sub-sample was 5 × 10⁻¹² erg s⁻¹ cm⁻² (in the 0.1 − 2.4 keV band), which excluded A750, A2187, A2396, A2631, and A2645 from the full HeCS sample. We then excluded three clusters from the flux-limited sample that were found to be clearly dominated by AGN upon inspection of available X-ray observations (these were A689 and A1366 based on Chandra data, and A2055 based on XMM-Newton data). More specifically, A1366 and A2055 have BLLacs in the brightest cluster galaxy (BCG), and A689 has a BLLac very close to the BCG (Giles et al. 2012). This resulted in a flux-limited sample of 50 clusters. We then obtained Chandra follow-up observations of those clusters lacking Chandra data.

To construct the final sample, we then excluded four clusters due to flaring in the available Chandra observations (A267, A667, A1361, and A1902), one cluster as it was part of a double cluster (A1758), and one cluster as it was off-chip in the Chandra observation available (A646 – the observation was of a radio galaxy in the cluster, but around 7′ from the cluster core). The final sample then consists of 44 clusters, which are summarised in Table 1. All clusters are in the redshift range 0.1 < z < 0.3. We call our sample the CHeCS (Chandra observations of Hectospec Cluster Survey clusters) sample. We note that the masses from Rines et al. (2013) were converted to be consistent with the value of H₀ we use in this paper.

3. Analysis

3.1. X-ray data

The Chandra data analysis follows closely that presented in Giles et al. (2015), which in turn is based on the analysis presented in Maughan et al. (2008, 2012). We summarise the main steps here, and highlight any differences in the updated analysis used in this paper. Any peculiarities of the analysis of individual clusters are noted in Appendix A.

All of the clusters in our sample were analysed with the CIAO⁴ 4.10 software package and CALDB⁵ version 4.8.1 (Fruscione et al. 2006). Projected emissivity profiles of the ICM were obtained from the surface brightness profile, measured in circular annuli centred on the centroid of the X-ray emission. The centroid location was measured in an X-ray image produced in the 0.7 − 2 keV band, and was determined by iteration using a circular region of 150″ and then 50″ initially centred on the X-ray peak (Maughan et al. 2008). Projected temperatures profiles were obtained by extracting spectra in circular annuli also centred on the X-ray centroid. The spectra were fit with an absorbed APEC plasma model (Smith et al. 2001) in XSPEC (Arnaud 1996). For the spectral analysis, instead of using the χ² statistic as in Maughan et al. (2016), we used the C-statistic, which assumes Poissonian statistics, and the spectra were grouped to contain at least five counts per energy bin. This minimal binning is required to avoid a bias in the XSPEC implementation of the C-statistic when there are bins in the spectrum with zero or very few counts (see Willis et al. 2005, and the discussion in Appendix B in the XSPEC manual). We exclude the outermost temperature measurement in the profile if the measured temperature in the final bin is greater than 15 keV, has errors > 50%, or including it leads to an unphysical increase in the best fit temperature profile.

A further minor difference from our previous analysis is that the absorbing column used in the spectral fits are fixed at the NHTOT value from Willingale et al. (2013), instead of the absorbing column from Dickey & Lockman (1990). The updated NHTOT value now contains absorption contributions from molecular hydrogen, not previously accounted for. We note this update does not significantly affect our results.

In order to calculate hydrostatic mass profiles, we use the models for the 3D temperature and gas density profiles as presented in Vikhlinin et al. (2006) and fit them to the observed projected temperature profile and surface brightness profile respectively (this mass fitting method is sometimes referred to as the forward-fitting mass method e.g. Ettori et al. 2013 and more details of our implementation are given in Appendix B). In this work we fit the 3D temperature and gas density models simultaneously, in contrast to the X-ray mass method used in Maughan et al. (2016). This fitting is done using a Markov chain Monte Carlo (MCMC) method using emcee (Foreman-Mackey et al. 2013), and the priors used on the 3D gas density and temperature models are presented in Table 2. Uncertainties on the mass profile are obtained straightforwardly from the MCMC chains. We also compute the mean R₅₀₀, and mean hydrostatic mass at R₅₀₀ from these chains (presented in Table 3). Gas mass profiles can also be obtained from the 3D gas density fit, with uncertainties obtained from the MCMC chains in the same way. This means that the uncertainties derived on the gas mass within R₅₀₀ self-consistently include the uncertainties on R₅₀₀.

To investigate the influence of the clusters’ dynamical states on our results, we split our sample into relaxed cool core (RCC) and non-relaxed cool core (NRCC) clusters, using the same three diagnostics as in Maughan et al. (2016). To be classed as RCC, a cluster has to have a low central cooling time (< 7.7 Gyr), a highly peaked central density profile (logarithmic slope > 0.7 within 0.048 R₅₀₀), and a low centroid shift (< 0.009). This definition of what constitutes a RCC cluster is a conservative one: for our sample of 44 clusters, 10 are classified as RCC (see Table 3).

This work is an extension of Maughan et al. (2016) from 16 to 44 clusters. As noted Sect. 2, we have dropped A267 and A2631 from the present analysis, leaving 14 clusters in common. For these clusters, we compared the hydrostatic masses at R₅₀₀ between the two analyses, finding a weighted average ratio of the new masses to those from Maughan et al. (2016) of 0.96 ± 0.05. This demonstrates that the improvements in our analysis, and changes to the Chandra calibration have not significantly affected our hydrostatic mass estimates.

3.2. Caustic masses

In this work, we use the caustic masses from HeCS, as presented in Rines et al. (2013). We summarise the key aspects of the analysis in the following. The caustic mass method uses spectroscopic redshifts for a number of member galaxies in a cluster (on average ∼180 for our sample, see Col. 7 Table 1), to determine their line-of-sight velocity (relative to the median of the distribution of the line-of-sight velocities of the galaxies of the main group on the binary tree Diaferio 1999), which when combined with their projected distances from the cluster centre, can be used to estimate a mass profile of a cluster.

By plotting the line-of-sight velocity versus the projected distance from the cluster centre for each of the member galaxies (referred to as a redshift diagram), an overpopulated region of this parameter space is clearly seen (see e.g. Rines et al. 2013), and is bounded by a ‘trumpet shape’. The boundaries between galaxies within this trumpet shape and outside it are known as caustics. The galaxies within this overpopulated parameter space are assumed to be bound by the cluster’s gravitational potential, and their velocities are assumed to be less than the escape velocity of the cluster. The amplitude, 𝒜(r), of this trumpet shape is related to the cluster’s escape velocity and the velocity anisotropy parameter β(r), and decreases as a function of projected distance from the cluster centre. Diaferio & Geller (1997) show that this caustic amplitude, 𝒜(r), can be related to the cluster mass within a radius, r:

$$\begin{array}{c}\hfill GM(<r)={\mathcal{F}}_{\beta }{\displaystyle {\int }_0^r{\mathcal{A}}^2(r)\mathrm{d}r.}\end{array}$$(1)

Equation (1) simplifies the correct equation $GM(<r)={\int }_0^r{\mathcal{F}}_{\beta }(r){\mathcal{A}}^2(r)\text{d}r$ by replacing the function ℱ_β(r) with the constant filling factor ℱ_β. The function ℱ_β(r) combines the cluster density profile and the anisotropy of the velocity field, both of which vary as a function of r. Beyond ∼0.5 R₂₀₀, ℱ_β(r) is roughly constant with radius, so assuming a constant filling factor ℱ_β should not lead to significant overestimation or underestimation of mass at larger radii (Diaferio 1999). However, within ∼0.5 R₂₀₀, approximating ℱ_β(r) to a constant might lead to an overestimate of the cluster mass by ∼15% (Serra et al. 2011).

We note that we use ℱ_β = 0.5 for the filling factor, which is appropriate for the algorithm of Diaferio (1999) used to calculate the caustic masses presented in Rines et al. (2013).

The uncertainties on the caustic masses represent the statistic precision of the caustic mass estimate due to the uncertainties in the location of the caustics. They do not represent the accuracy with which the caustic mass measures the true mass, which is ≈30% (Serra et al. 2011), driven largely by the effects of projection and the viewing angle of the cluster (this is accounted for in our model below by an intrinsic scatter between caustic and true masses).

The uncertainty on the caustic location increases with decreasing ratio ρ between the number of galaxies within the caustics and the number of galaxies in the redshift diagram. On average, this recipe is a robust estimate of the statistical uncertainty (Serra et al. 2011). However, in sparsely-sampled redshift diagrams, the caustic technique might underestimate the statistical uncertainty, because the ratio ρ might remain relatively large, despite the low number of galaxies in the diagram. We verified that these cases have a negligible effect on our results by repeating our analysis with a uniform 50% error on all caustic masses. In this test, the intrinsic scatter term in our model is effectively removed and we prevent any points with underestimated errors from driving the results.

The number of galaxies within the caustics, Ngal, can also have a systematic effect on the mass derived using the caustic method. The simulations of Serra et al. (2011) suggest an underestimate of the mass of ∼10% at R₅₀₀ for Ngal ∼ 100, with the masses being unbiased for 200 ≲ Ngal ≲ 1000. Any underestimate of mass due to low values of Ngal is expected to approximately cancel with the overestimate of mass at R₅₀₀ due to the assumption of constant filling factor. However, as discussed below (Sect. 4.2), our results suggest that the magnitude of the mass underestimate for low values of Ngal may be larger than predicted, and our main conclusions are thus based on the clusters with the largest values of Ngal.

3.3. Modelling the mass biases

In order to constrain the bias and scatter between the hydrostatic and caustic mass measurement techniques, we use a Bayesian framework (of which a graphical representation is shown in Fig. 1). The modelling framework is the same as in Maughan et al. (2016), but we summarise the key points here. In the following, we use μ to denote the base 10 logarithm of mass.

Fig. 1. Graphical model summarising the mass modelling framework that we present in Sect. 3.3. The symbols are also described in that section; in summary, the symbol μ represents mass, κ is the bias term, δ is the scatter term and σ is the measurement error term (all in log space), and the C and X subscripts denote whether the quantity refers to that from the hydrostatic or caustic method. The refers to an observed quantity. 𝒩 refers to a normal distribution, and ‘∼’ denotes ‘is distributed as’.

A given cluster with a ‘true’ mass, μ also has a ‘true’ caustic mass, μ_C, and ‘true’ hydrostatic mass, μ_X, which are related to the true mass, μ, as follows:

$$\begin{array}{c}\hfill {\mu }_X\sim \mathcal{N}(\mu +{\kappa }_X,{\delta }_X)\end{array}$$(2)

$$\begin{array}{c}\hfill {\mu }_C\sim \mathcal{N}(\mu +{\kappa }_C,{\delta }_C)\end{array}$$(3)

where κ_X and κ_C parametrise the bias between the true mass and the ‘true’ hydrostatic and caustic mass, respectively, and δ_X and δ_C parametrise the intrinsic scatter between the true mass and the ‘true’ hydrostatic and caustic mass, respectively. 𝒩 refers to a normal distribution, and “∼” is equivalent to “is distributed as”.

These ‘true’ hydrostatic and caustic masses are related to the observed hydrostatic and caustic masses ($\widehat{{\mu }_X}$ and $\widehat{{\mu }_C}$, respectively), as follows:

$$\begin{array}{c}\hfill \widehat{{\mu }_X}\sim \mathcal{N}({\mu }_X,{\sigma }_X)\end{array}$$(4)

$$\begin{array}{c}\hfill \widehat{{\mu }_C}\sim \mathcal{N}({\mu }_C,{\sigma }_C)\end{array}$$(5)

where σ_X and σ_C represent the standard deviation of the lognormal error for the observed hydrostatic and caustic masses, respectively. Note that our model does not take into account selection effects that could arise due to the X-ray flux limit. However, we show in Sect. 5.3 that this has a negligible impact on our results.

We used the same (loose) priors as in Maughan et al. (2016) on the cluster masses, bias and scatter terms. Specifically, for the cluster masses in log space, μ, we applied a uniform prior from 12 to 17. For the bias terms in log space, κ_X and κ_C, we applied normal priors with mean 0 and a standard deviation of 1. For the intrinsic scatter terms in log space, δ_X and δ_C, we applied normal priors (truncated at zero) with mean 0.09 and a standard deviation of 2.2, corresponding to a weak prior with a mean of ∼20% in normal space.

Using the model described above, we used the observational data that we have for each cluster ($\widehat{{\mu }_X}$, σ_X, $\widehat{{\mu }_C}$, σ_C) to constrain the ‘true’ hydrostatic and caustic mass terms (μ_X, μ_C), the scatter terms (δ_X and δ_C) and bias terms (κ_X and κ_C) for all of the clusters in our sample. Both the scatter terms and the bias terms are degenerate, however, the intrinsic scatter between the hydrostatic and caustic mass measurements

$$\begin{array}{c}\hfill \delta =\sqrt{{\delta }_X^2+{\delta }_C^2}\end{array}$$(6)

and the mean bias between the hydrostatic and caustic mass measurements

$$\begin{array}{c}\hfill \kappa ={\kappa }_X-{\kappa }_C={\mu }_X-{\mu }_C={log}_{10}\left(\frac{M_X}{M_C}\right)\end{array}$$(7)

can be constrained by the data.

When reporting the values of κ and δ in Tables 4 and 5, we report the median value of δ, with errors given as the difference between the median and the 16th and 84th percentiles. For κ, we report the mean value. As κ = log₁₀(M_X/M_C), and κ (a quantity in log space) is normally distributed, the posterior of M_X/M_C (a quantity in linear space) is lognormally distributed. Therefore, in Tables 4 and 5 we summarise the posterior of M_X/M_C by reporting its median and errors in the same way that we did for δ. In these tables, the κ values are given in base 10 log space, and the δ values have been converted to percentage scatters.

We also note that in Table 3 and in all figures showing mass errors, we report the mean and standard deviation of the hydrostatic and caustic mass measurements in linear space, which is the form in which masses are often reported in other work. However, we find that the uncertainties on the mass measurements are in fact better described using lognormal errors, and so choose to model the biases associated with the caustic and X-ray hydrostatic mass measurement methods in log space.

We use the probabilistic programming language STAN⁶ to implement the model described above, specifically using the No-U-Turn Sampler (Hoffman & Gelman 2011). We sampled the parameters in our model with four chains of 5000 steps each. The analysis can be done at any radius, and so by repeating the analysis at increasing radii, we produced a profile of the mean bias between the two mass measurement methods.

In our analysis, the mass determination methods are independent and use different definitions of the cluster centre (the X-ray centroid for the hydrostatic mass, and the peak of the galaxy distribution on the sky and the median of the line-of-sight velocity distribution of the cluster members with the caustic method). This means that the radial coordinates are not identical and a given radius does not correspond to an identical aperture for both methods. For our main analysis, we scale radial coordinates to the R₅₀₀ estimated from the hydrostatic method. This is motivated by the expectation that the hydrostatic mass should have lower intrinsic scatter with the true mass. We show in Sect. 5.2 that using R₅₀₀ computed independently for each method does not alter our conclusions.

4. Results

In the following, we first present the comparisons between hydrostatic and caustic masses for the full CHeCS sample. We then demonstrate that the caustic masses for clusters with lower values of Ngal are likely to be underestimated. Finally we present our main results, based on the subset of clusters with higher values of Ngal.

4.1. Initial results for the full sample

The caustic and hydrostatic mass profiles for each cluster are shown in Fig. C.1⁷. These profiles were used to calculate the M_500, X and M_500, C, using the R₅₀₀ value calculated from the hydrostatic mass profile for each cluster (unless stated otherwise, R₅₀₀ always refers to that derived from the hydrostatic mass profile)⁸. The hydrostatic and caustic masses are reported in Table 3 and are compared in Fig. 2.

Fig. 2. The hydrostatic versus the caustic masses are plotted for the full sample of clusters (with masses calculated at the hydrostatic R500). NRCC clusters are red circles and RCC clusters are blue triangles. The solid black line shows the best fitting model to the full data, with the shaded region indicating the 1σ uncertainty. The 1:1 line is also plotted as a dashed black line.

The temperature profiles required for the hydrostatic mass estimates were measured out to radii close to, or beyond, R₅₀₀ for the majority of clusters (see Fig. C.1). For the mass profiles at radii greater than the extent of the temperature profile, we extrapolated the best fitting 3D temperature profile. As a consequence, our hydrostatic mass profiles beyond the extent of the measured temperature profiles are less robust. The median radius out to which the temperature profiles were measured is 0.95 R₅₀₀, and the range of radii is 0.51−1.62 R₅₀₀.

The use of the hydrostatic estimate of R₅₀₀ to determine both sets of masses introduces a covariance between the hydrostatic and caustic masses that is not included in our model. For this reason, we verified that using a fixed aperture of radius 1 Mpc for the mass measurements made a negligible difference to our results (as found in Maughan et al. 2016).

The M_X/M_C profile of each cluster (calculated as $\widehat{{\mu }_X}-\widehat{{\mu }_C}$) is shown in Fig. 3. We also plot the mean M_X/M_C profile (calculated as the mean bias κ, see Eq. (7)). The M_X/M_C profile is consistent with the bias being constant with radius. At smaller radii the M_X/M_C value decreases, but this is likely due to the caustic masses being overestimated at small radii (see e.g. Fig. 12 Serra et al. 2011).

Fig. 3. The MX/MC profiles for all 44 clusters in our sample, scaled in radius by the hydrostatic R500. NRCC clusters are solid red lines and RCC clusters are blue dot dashed lines. The 1:1 line is also plotted as a thin dashed black line. The thick dashed black line shows the average MX/MC ratio, and the shaded region shows the 1σ uncertainty.

The M_500, X/M_500, C ratios at R₅₀₀ for the individual clusters is shown in Fig. 4, and the average ratios are summarised in Table 4. At R₅₀₀ the caustic and hydrostatic masses generally do not agree well. For the full sample, the hydrostatic mass on average ∼30% higher than the caustic mass (significant at ≈3σ). Similar results are found for both the NRCC and RCC clusters.

Fig. 4. We show the M500, X/M500, C ratio for all 44 clusters in our sample. NRCC clusters are red circles and RCC clusters are blue triangles. The solid black line and shaded region show the average M500, X/M500, C ratio at the hydrostatic R500. 1σ errors are shown. The clusters are ranked by Ngal with the highest at the top, and the horizontal lines separate the clusters into the low, mid and high Ngal subsets (see Sect. 4.2).

This M_500, X/M_500, C value is consistent with the our earlier findings (Maughan et al. 2016), but rules out M_500, X/M_500, C = 1 at higher significance, and remains counter-intuitive. The hydrostatic mass is expected to be underestimated compared to the true mass (e.g. Nelson et al. 2014; Eckert et al. 2019) while the caustic mass expected to be overestimated compared to the true mass at R₅₀₀ (e.g. Serra et al. 2011); thus we expect the M_500, X/M_500, C ratio to be less than 1.

4.2. Dependency on Ngal

Figure 2 suggests that the high value of the M_500, X/M_500, C ratio for the full sample may be driven by clusters with lower caustic masses (with hydrostatic masses being overestimated and/or caustic masses being underestimated). However, simulations indicate that dynamical mass estimators in general (including a different implementation of the caustic method to that used here) tend to underestimate masses when the number of galaxies used is low (Wojtak et al. 2018). This is a confounding factor since the lower mass clusters tend to have lower Ngal.

In order to investigate the role of Ngal in our results, the sample was split on Ngal into three approximately equal subsets of 15, 15 and 14 clusters (see the Ngal sub-sample column in Table 3). These low, mid and high Ngal sub-samples have 63 ≤ Ngal ≤ 129 (mean 93), 131 ≤ Ngal ≤ 208 (mean 174), and 210 ≤ Ngal ≤ 461 (mean 281) respectively. Figure 5 shows the comparison between hydrostatic and caustic masses with the different Ngal subsets indicated. This illustrates the trend that lower mass clusters tend to have lower Ngal. It is also clear that, at a given hydrostatic mass, the clusters with the lowest Ngal are most discrepant in their caustic mass.

Fig. 5. The hydrostatic and caustic masses measured within the hydrostatic R500 are plotted for the full sample. This is the same data as plotted in Fig. 2, but now showing the Ngal subsets. The high, mid and low Ngal bins are green squares, blue triangles and red circles respectively. The 1:1 line is also plotted as a dashed black line.

To investigate this further, the gas fraction, f_gas (defined as the gas mass, m_gas, within some radius divided by the total mass within the same radius), was found to be a useful discriminator between possible systematics in the analysis. This is because m_gas can be measured accurately from X-ray observations (e.g. Nagai et al. 2007), and there are robust theoretical predictions for the range of f_gas expected. We thus computed f_gas at a radius of 1 Mpc (chosen to avoid dependency of the radius on any mass estimation method) with either the hydrostatic and caustic masses in the denominator. The results are plotted in Fig. 6. The plot also indicates the approximate range of f_gas (∼0.1−0.15) expected at this radius (e.g. Eckert et al. 2016).

Fig. 6. Left: fgas measured at 1 Mpc using caustic mass in the denominator is plotted against Ngal. Right: same plot, but with fgas calculated with hydrostatic mass in the denominator. In both panels, the black dashed lines indicate fgas = 0.10 and fgas = 0.15 for reference, encompassing the approximate expected range of fgas.

When the caustic mass is used as the denominator in f_gas, there is a strong indication of a trend with Ngal, with unphysically high values of f_gas for the clusters with lower values of Ngal. The mean and standard error of f_gas was 0.23 ± 0.03, 0.15 ± 0.01 and 0.15 ± 0.01 for the low, mid and high Ngal sub-samples respectively.

In contrast, when the hydrostatic mass is used to compute f_gas, the trend with Ngal disappears, and the average values of 0.12 ± 0.01, 0.11 ± 0.01, 0.12 ± 0.01, for the low, mid and high Ngal sub-samples respectively are consistent with the expected range.

These f_gas measurements suggest that Ngal is responsible for the high values of M_500, X/M_500, C due to an underestimate of M_C. However, due to the correlation between Ngal and mass, the possibility remains that a trend with mass is more fundamental. During the completion of the work presented here, additional caustic mass measurements became available in the HeCS-Omnibus catalogue (Sohn et al. 2020) that allow the effects of mass and Ngal to be separated. The new catalogue includes new caustic mass estimates for clusters in our sample, using the same caustic method, but based (in most cases) on a larger number of galaxies per cluster.

When comparing the earlier and newer caustic mass measurements, it was found that for well-sampled clusters, increasing the sampling had little effect on the mass. On the contrary, for more poorly sampled clusters, increasing Ngal systematically led to an increased mass estimate. This is illustrated in Fig. 7, which shows the ratio of the caustic masses (measured at the hydrostatic R₅₀₀) from Sohn et al. (2020) to those from both Rines et al. (2013) (as used in our work), plotted against the difference in Ngal between the two analyses.

Fig. 7. The ratio of the caustic masses (measured at the hydrostatic R500) from Sohn et al. (2020) to those from Rines et al. (2013) are plotted as a function of the difference in Ngal between the analyses. The clusters are split into subsets containing those where the minimum Ngal used in the two analyses was at least 200 (green squares) and those where the minimum Ngal was less then 120 (red circles). The lines show the best fits to each subset.

The clusters were divided into subsets based on the minimum Ngal used in the two caustic calculations. For clusters where at least 200 galaxies were used in both analyses (green squares in Fig. 7), increasing the sampling has a weak or null effect on the mass estimate. The green line is a straight line fit to the data in this parameter space (i.e. log(M_Sohn/M_Rines) versus N_Sohn − N_Rines) using BCES y on x regression (Akritas & Bershady 1996). This is a simplistic model and is not intended to fully describe the systematic effects in the analyses, but is an adequate description of the data presented here. The slope of the line is 0.001 ± 0.001.

The impact of increasing Ngal for less well-sampled clusters is seen in the blue triangles in Fig. 7. These are clusters for which Ngal is < 120 in at least one of the two analyses. It is clear in the figure that increasing the sampling leads to systematically higher caustic mass estimates. The best fitting line to this subset has a slope of 0.008 ± 0.002. The subsets were defined to split the sample into three equal groups, and the intermediate subset falls between the two outer cases and are omitted from the figure for clarity.

We further test the role of Ngal as follows. We undersample the galaxy catalogues of the 14 clusters in our high-Ngal sub-sample, by progressively removing the faintest galaxies. When Ngal drops by ∼50%, falling in the range ∼100 − 200, the caustic mass within R₅₀₀ decreases by ∼15%, on average. This average mass underestimate is consistent with the analysis of Serra et al. (2011) who found that with Ngal ∼ 100 the escape velocity at R₅₀₀ is underestimated by ∼5 − 10%, whereas Ngal ≳ 200 is required to return unbiased escape velocity estimates (see also Halenka et al. 2022). The lower Ngal of the undersampled clusters in our high-Ngal sample would move these clusters in our mid-Ngal sample: the ∼15% mass underestimate we find in this test is thus also consistent with the larger ratio M_500, X/M_500, C we find for our mid-Ngal sample, as illustrated in Sect. 4.3 below (Table 5).

Based on these results, we conclude that a systematic underestimate of the caustic mass at R₅₀₀ for low values of Ngal is responsible for the high apparent M_500, X/M_500, C ratio found at low masses, and that the high Ngal sub-sample (Ngal ≥ 210) is not significantly affected by this systematic. We thus repeat our analysis on the separate Ngal subsets in the following section, and will draw our final conclusions based on the results for the high Ngal subset alone.

4.3. Final results

The analysis described in Sect. 4.1 was repeated for each of the Ngal subsets separately, and the results are summarised in Table 5. Figure 8 shows the profiles of M_X/M_C for the different Ngal subsets. The separation between the subsets is clear, with no radial dependency.

Fig. 8. The average MX/MC profiles for the high Ngal, mid Ngal and low Ngal clusters in green lines, blue lines, and red dot-dashed lines respectively. The 1:1 line is also plotted as a thin black dashed line. The thick dashed black line shows the average MX/MC profile for all clusters. The profiles are scaled in radius by the hydrostatic R500, and the shaded regions show the 1σ uncertainties.

As expected based on the preceding analysis, the results in Table 5 show a trend towards lower M_X/M_C ratios with increasing Ngal. The most reliable results are those based on the sample of 14 clusters with the highest values of Ngal, and these are highlighted in boldface in the table. For this subset, the average ratio of hydrostatic to caustic masses is $M_{500,X}/M_{500,C}=1.{12}_{-0.10}^{+0.11}$ with an intrinsic scatter between the masses of $\delta =({25}_{-08}^{+12})\%$.

The comparison between the hydrostatic and caustic masses for the high Ngal subset is shown in Fig. 9, with the RCC/NRCC classification of the clusters indicated. There does not appear to be a significant trend with dynamical state but given the reduced size of this high Ngal subset, a more detailed analysis would not be useful.

Fig. 9. We show the hydrostatic masses versus the caustic masses for the 14 clusters in our sample with the highest Ngal values, both calculated at the hydrostatic R500. We show RCC clusters as blue triangles and NRCC clusters as red circles. The solid black line shows the best fitting model, with the shaded region indicating the 1σ uncertainty. The 1:1 line is also plotted as a dashed black line.

Recall that 14 of the 16 clusters used in our earlier analysis in Maughan et al. (2016) are in common with the current analysis. Of these, six are in our high Ngal subset, seven are in the mid Ngal sub-sample, and one is in the low Ngal sub-sample. This mix of mid- and high-Ngal clusters naturally explains why Maughan et al. (2016) found an M_X/M_C ratio that was slightly higher (M_500, X/M_500, C = 1.20 ± 0.12) than our final result for the high Ngal subset.

5. Discussion

In the following section we consider the possible systematics that could have impacted our results. We explore in more detail our conclusion that the apparent increase in M_500, X/M_500, C at low masses is driven by Ngal rather than by mass (or some other factor). Finally, we compare our results with other related analyses, and discuss the constraints that can be derived on the level of hydrostatic bias in massive clusters.

5.1. The impact of Ngal

Our investigation of the comparison between hydrostatic and caustic masses suggested that the caustic masses were underestimated when fewer galaxies were used. While a general result of dynamical masses being biased low when based on smaller numbers of galaxies was shown in Wojtak et al. (2018), the caustic method examined in that work differed from the implementation used here.

The dependence on Ngal of the caustic mass was investigated for our implementation of the caustic method by Serra et al. (2011). They found that the true escape velocity profile will be underestimated by ∼5−10% at R₅₀₀ (Fig. 20, Serra et al. 2011) for clusters with Ngal ∼ 100 (compared to clusters with a higher Ngal for which the true escape velocity is not underestimated). This would lead to the caustic mass being underestimated by ∼10−20% at R₅₀₀, were it not for the assumption of a constant ℱ_β (see Sect. 3.2), which leads to a ∼10−15% overestimate of the true mass at R₅₀₀ for clusters with a higher Ngal (Fig. 12, Serra et al. 2011). Thus, the simulations suggest that these two effects approximately cancel out at our comparison radius of R₅₀₀, such that a large underestimation of the cluster mass using the caustic method for low Ngal clusters (like that implied by our results) is not expected.

As discussed in Sect. 4.2, we consider the caustic masses for the high Ngal subset to be unaffected by this apparent bias, and base our main conclusions on this subset. As a test of the robustness of our results, we repeated the measurement of M_500, X/M_500, C using the HeCS-Omnibus caustic data of Sohn et al. (2020) in place of the original HeCS caustic data of Rines et al. (2013) that was used for our main analysis. We selected the 14 clusters in our sample with the highest Ngal in HeCS-Omnibus (Ngal > 274 with a median of 413). For these best-sampled clusters, using the HeCS-Omnibus caustic masses gives an average hydrostatic to caustic mass ratio of M_500, X/M_500, C = 1.10 at the hydrostatic R₅₀₀. This agrees very well with the M_500, X/M_500, C = 1.12 we found for the 14 clusters in our high Ngal subset when using the Rines et al. (2013) data (for which Ngal > 210 with a median of 234). This supports our conclusion that our high Ngal subset is unaffected by systematics due to low numbers of galaxies (i.e. Ngal ≳ 200 is sufficient to avoid biases in the caustic mass at R₅₀₀). Indeed, when we apply the caustic technique to the cluster catalogues in the high-Ngal sample where we undersampled Ngal such that Ngal ≲ 200, the caustic mass within R₅₀₀ is underestimated by ∼15%, on average, suggesting that the increase of the ratio M_500, X/M_500, C from ∼1.12 of the high-Ngal sample to ∼1.28 of the mid-Ngal sample is mostly driven by the lower number of the caustic members.

Poorer sampling of the caustics is likely to be responsible (at least in part) for some of the outliers seen in the comparisons between the hydrostatic and caustic masses (e.g in Fig. 4). When restricted to the high Ngal subset, there are no significant outliers in the comparison between hydrostatic and caustic masses, given the ∼30% intrinsic scatter between the mass estimators (e.g. Fig. 9).

5.2. Cluster centres and apertures

In the caustic analysis, the cluster centre is located as part of the process (Rines et al. 2013), while the centre is fixed at the X-ray centroid for the hydrostatic analysis. This means that the profiles are not generally centred at the same location. The separation between the centres used for the X-ray and caustic analyses were generally small, and showed no trend with Ngal or with M_500, X/M_500, C (the median separation was ≈130 kpc in each of the Ngal subsets). We thus conclude that the increasing disagreement between hydrostatic and caustic masses at low Ngal is not being driven by miscentring.

Serra et al. (2011) investigated the impact of miscentering of the caustic analyses on the mass determination in their simulations. At R₅₀₀, the caustic mass estimate was increased by ≈10% when the profile was forced to use the true centre. This suggests that our measurements of M_500, X/M_500, C could be overestimated by ≈0.1 due to miscentring of the profiles.

When the separation between the centres used for the caustic and hydrostatic analyses is small, then the use of a common aperture size (i.e. the hydrostatic R₅₀₀) allows the mass estimators to be compared without adding scatter from different aperture definitions. However, given that the centres are not identical, it is also helpful to compare the masses when they are derived fully independently (i.e. comparing the caustic mass within the caustic R₅₀₀ with the hydrostatic mass within the hydrostatic R₅₀₀).

We repeated our analysis using these independent definitions of R₅₀₀. For the full sample of 44 clusters, we found that mass ratio increased from $M_{500,X}/M_{500,C}=1.{33}_{-0.09}^{+0.10}$ (using the hydrostatic R₅₀₀) to $M_{500,X}/M_{500,C}=1.{52}_{-0.14}^{+0.16}$ (using independent R₅₀₀). In the case of the high Ngal sub-sample, the original $M_{500,X}/M_{500,C}=1.{12}_{-0.10}^{+0.11}$ (hydrostatic R₅₀₀) increased slightly to $M_{500,X}/M_{500,C}=1.{16}_{-0.13}^{+0.16}$ (independent R₅₀₀).

These results make sense. When the masses agree well (as in the high Ngal sub-sample) then the R₅₀₀ definitions also agree so the results do not change. If the masses are more discrepant (as in the full sample) then using the independent R₅₀₀ will tend to amplify the difference.

As might be expected, the scatter between the mass estimators is somewhat larger when the independent R₅₀₀ are used. The scatter increases from (39 ± 6)% to $({55}_{-8}^{+9})\%$ for the full sample, and from $({25}_{-08}^{+12})\%$ to $({35}_{-11}^{+15})\%$ for the high Ngal sub-sample.

Based on this analysis, we conclude that the use of the hydrostatic R₅₀₀ for both mass determinations does not have a strong effect on our results.

5.3. Impact of selection effects

In our work we do not account for the impact of selection effects in our model. The sample is selected on the basis of X-ray flux, so covariance between X-ray flux and M_X could have an impact on the inferred M_500, X/M_500, C ratio. Specifically, if there were a strong positive covariance between X-ray luminosity, L_X, and hydrostatic mass, M_X, given the true mass, then X-ray selection would select clusters that were both more luminous than average and had higher hydrostatic masses than average, given their true mass. The size of this effect would depend both on the scatter in the relevant relations and the degree of covariance between M_X and L_X.

We assessed the possible impact of this selection effect by simulating cluster populations with different degrees of covariance between L_X and M_X, imposing the sample selection function, and then applying our model (described in Sect. 3.3) to the simulated population to determine the impact of the sample selection on the recovered M_500, X/M_500, C ratio.

In more detail, we sampled a large number of points from a mass function where each point consists of a true mass, M_true, and a redshift. We use the mass – luminosity relation from Mantz et al. (2010) to assign the true luminosity, L_true to each cluster. We then assume a scatter between L_true and L_X of 30%, a scatter between M_true and M_X of 10% and a scatter between M_true and M_C of 30%. We also require a correlation coefficient, ρ_LM, between L_X and M_X to create the covariance matrix for the multi-dimensional Gaussian from which we sample L_X and M_X from L_true and M_true. It is the strength of this correlation coefficient that dictates how much a selection on X-ray flux might bias an inference based on hydrostatic masses, and we investigate the impact of different assumed values below. Finally, we apply our X-ray selection to the fluxes derived from the L_true and redshift values, and use our model to measure the M_500, X/M_500, C ratio for the resulting simulated sample. We do not set any hydrostatic or caustic bias in the simulation, so in the absence of selection biases, we would measure M_500, X/M_500, C = 1.

To our knowledge, covariances between L_X and hydrostatic mass have not been measured in the literature, but the covariances between other ICM properties have been investigated in a few studies (e.g. Maughan 2014; Mantz et al. 2016; Andreon et al. 2017b; Farahi et al. 2018). For example, for the covariance between L_X and m_gas, Mantz et al. (2016) found a correlation coefficient of 0.43, while Farahi et al. (2018) found 0.76. Meanwhile, for the covariance between L_X and T, Mantz et al. (2016) found a correlation coefficient of 0.53, while Farahi et al. (2018) found 0.49. In both cases, the values obtained by Mantz et al. (2016) are the more relevant here since they make use of a soft-band, core-included X-ray luminosity, which matches the selection method used for our sample (Farahi et al. 2018, made use of bolometric luminosities with the core regions removed, but we include their results here for comparison). We therefore conclude that a correlation coefficient of ρ_LM = 0.5 is a reasonable upper limit on the strength of the covariance between L_X and M_X and use this as a reference value for our simulations. We also performed the simulations with with ρ_LM = 1 as a maximally pessimistic case.

We find that for ρ_LM = 0.5, the recovered M_500, X/M_500, C value is 1.03 ± 0.01 and for ρ_LM = 1 we measured M_500, X/M_500, C = 1.07 ± 0.01 (for a sample of 500 simulated clusters).

We also checked for a mass dependence of the selection bias by splitting the simulated clusters into low (log(M_X) < 14.5), mid (14.5 < log(M_X) < 14.75), and high (log(M_X) > 14.75) mass sub-samples (each containing 500 simulated clusters). The measured M_500, X/M_500, C ratio was largest for the low-mass sub-sample, due to a larger proportion of the clusters being close to the flux limit. For the low-mass sub-sample, assuming ρ_LM = 0.5 we found M_X/M_X = 1.06 ± 0.01, increasing to 1.15 ± 0.01 for ρ_LM = 1. For the high mass sub-sample (which is similar in mass to our high Ngal subset), the measured M_500, X/M_500, C ratios dropped to 1.02 ± 0.01 and 1.03 ± 0.01 for ρ_LM = 0.5 and ρ_LM = 1 respectively.

We conclude that, while neglecting selection effects tends to bias the M_500, X/M_500, C ratio high, the impact is not significant for the mass range of the high Ngal subset, and that selection biases alone cannot explain the high values of M_500, X/M_500, C that we find in the lower mass, lower Ngal clusters.

5.4. Impact of X-ray temperature calibration

It is established that there are systematic differences in cluster temperatures measured by Chandra and XMM-Newton, temperatures measured with Chandra being systematically hotter (Schellenberger et al. 2015). This has been found to have a corresponding systematic impact on the derived hydrostatic masses (Mahdavi et al. 2013; Schellenberger et al. 2015). We note, however, that Martino et al. (2014) found no systematic differences in hydrostatic masses for a sample of 50 clusters with both XMM-Newton and Chandra data.

Schellenberger et al. (2015) found that the hydrostatic masses measured from Chandra data were 14 ± 2% higher than those from XMM-Newton. They found no significant trend with cluster mass, suggesting that the apparent trend of M_500, X/M_500, C with mass or Ngal that we find is not driven by this systematic. If our hydrostatic masses were systematically overestimated by ≈15%, this would have a direct effect of reducing our derived M_500, X/M_500, C ratios by 0.15.

5.5. Summary of systematics

We have identified a number of possible systematics that could impact our comparison between caustic and hydrostatic masses. In the case of the caustic masses, poorer sampling of the caustics was found to result in an underestimate of the caustic masses for Ngal ≲ 200, producing an apparent trend of M_500, X/M_500, C with mass. For the better-sampled clusters, this systematic appears to be absent, but the assumption of constant ℱ_β is expected to lead to an overestimate of the caustic mass by 10 − 15% at R₅₀₀ (i.e. biasing M_500, X/M_500, C low by up to 0.15). The miscentering of the caustic profiles leads to an estimated underestimate of the caustic masses of ≈10% at R₅₀₀ (biasing M_500, X/M_500, C high by ≈0.1).

We find that the use of X-ray flux to select the sample combined with a possible covariance between X-ray luminosity and hydrostatic mass could bias the hydrostatic masses to be high. However, the effect is estimated at a few percent for the more massive clusters comprising our high Ngal subset, and the trend with mass is not sufficient to explain the strength of the trend of M_500, X/M_500, C with mass or Ngal. Finally, systematics on the calibration of Chandra could lead to the hydrostatic masses being overestimated by ≈15% (biasing M_500, X/M_500, C high by about 0.15).

5.6. Comparison with other work

Direct comparisons between hydrostatic and caustic masses for samples of clusters are scarce in the literature. Andreon et al. (2017a) compared the mean f_gas for a sample of 34 clusters with caustic masses, with the mean value from separate samples with hydrostatic masses from Chandra data. They inferred that caustic masses were slightly higher than, but consistent with, hydrostatic masses at R₅₀₀. This is consistent with the result we find for the high Ngal subset. However, the clusters with caustic masses used in Andreon et al. (2017a) cover a range in mass and Ngal that is comparable with our low and mid Ngal subsets, and are in tension (at the ≈2σ level) with those subsets for which we find the caustic masses to be lower than hydrostatic masses (although we consider our results for those subsets to be less reliable due to the poorer sampling of the caustics).

Foëx et al. (2017) compared hydrostatic and caustic masses for a sample of 13 clusters with masses, redshifts and Ngal similar to our high Ngal subset. They found a mean M_C/M_X = 1.32 ± 0.18 at R₅₀₀ (we note that this is the inverse of our mass ratio), and hence caustic masses were higher than hydrostatic masses by ∼30% (although the ratio dropped to unity when some substructures were removed in the dynamical analysis). While this contrasts with our results, the difference is not strongly significant (≈2σ). Furthermore, the Foëx et al. (2017) analysis was based on hydrostatic masses from XMM-Newton; adjusting these masses to match the Chandra calibration would reduce the tension with our result to less than 1σ.

Ettori et al. (2019) compared hydrostatic masses (from XMM-Newton) and caustic masses at R₂₀₀ for 6 clusters, finding that the caustic masses were lower than the hydrostatic masses by about 30% on average. This is in the same sense as the difference we find for our sample, but the clusters include some with lower Ngal than our high Ngal subset, so a direct comparison is not possible.

5.7. Implications for the hydrostatic bias

Conventionally, the hydrostatic bias is defined as b = 1 − M_X/M, where M_X is the hydrostatic mass and M is the true mass. Our best estimate of the ratio of hydrostatic to caustic masses, with statistical error is $M_{500,X}/M_{500,C}=1.{12}_{-0.10}^{+0.11}$, for the high Ngal subset. Under the strong assumption that the caustic masses are unbiased, we can write b = 1 − M_500, X/M_500, C, in which case our results exclude a hydrostatic bias of more than 20% at the 3σ level based on statistical uncertainty alone. However, taking the possible systematic uncertainties on both the caustic and hydrostatic masses into account, we cannot rule out larger values of the hydrostatic bias.

Constraints on the hydrostatic bias from direct comparisons between hydrostatic and WL masses currently encompass a range from approximately zero (e.g. Gruen et al. 2014; Israel et al. 2014; Applegate et al. 2016; Smith et al. 2016) up to around 30% (e.g von der Linden et al. 2014; Donahue et al. 2014; Hoekstra et al. 2015). Meanwhile, constraints on the level of bias have also been estimated indirectly by comparisons of WL masses with masses derived from scaling relations of other observables that were calibrated to hydrostatic masses (e.g. from X-ray luminosity Sereno & Ettori 2017 or the Sunyaev-Zel’dovich effect signals Medezinski et al. 2018; Miyatake et al. 2019). These also tend to favour biases of 20 − 30%.

The question of hydrostatic bias has also been addressed with simulations, which typically find biases of 10 − 30%, but with radial and mass dependences, and sensitivity to the method used to reconstruct the hydrostatic masses from the simulation data (e.g. Nagai et al. 2007; Rasia et al. 2012; Pearce et al. 2020; Barnes et al. 2021).

Simulations suggest a significant mass dependence of the hydrostatic bias, with the bias increasing for higher mass clusters (e.g Barnes et al. 2021). However, due to the relatively poorer sampling of the caustics in our lower mass clusters, we are not able to test these results observationally.

Overall, there is not a clear consensus on the magnitude of the hydrostatic bias. Our results are more consistent with those estimates from the low end of the range and add a useful new datapoint to this still-open question, since the caustic mass determination is subject to different systematics and assumptions to those affecting the WL masses which are primarily used to address this question observationally.

6. Summary and conclusions

We compared masses of 44 galaxy clusters determined independently from X-ray data under the assumption of hydrostatic equilibrium, and galaxy velocities using the caustic method. Our initial results showed significant discrepancies between the masses which appeared to increase for lower mass clusters. We found that this was driven by the poorer sampling of the caustics for the lower mass clusters leading to an underestimate of the mass using the caustic method at R₅₀₀. The caustic method is not optimised for mass determinations at these radii, but this finding suggests that simulations have underestimated the impact of poorer sampling at smaller radii.

Restricting our analysis to the 14 clusters with best-sampled caustics (the high Ngal subset), we find an average ratio of the hydrostatic to caustic masses of $M_{500,X}/M_{500,C}=1.{12}_{-0.10}^{+0.11}$. We find no evidence for radial dependence of this ratio or dependence on dynamical state of the cluster, but these conclusions are limited by the reduced sample size.

We investigated possible systematics that could affect the mass estimates for this high Ngal subset and find that the most significant are as follows. The assumption of constant ℱ_β and possible miscentering of the caustic profiles could bias the caustic masses high and low respectively by 10 − 15%. The uncertainties on the absolute calibration of the X-ray temperature measurements could bias the hydrostatic masses high by about 15% compared to those measured with XMM-Newton.

We interpret our results as favouring a hydrostatic bias that is close to zero. However, the systematics we identified would allow for significantly larger levels of hydrostatic bias.

Acknowledgments

C.L., B.J.M. and R.T.D. acknowledge support from the UK Science and Technology Facilities Council (STFC) grant ST/R00700/1. B.J.M. acknowledges further support from STFC grant ST/V000454/1. C.L. also acknowledges support from an ESA Research Fellowship. A.D. acknowledges partial support from the Italian Ministry of Education, University and Research (MIUR) under the Departments of Excellence grant L.232/2016, and from the INFN grant InDark. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: Notes on individual clusters

In this section we detail where we differed from the analysis methods presented in Sect. 3, or any particular points of interest for each cluster in our sample.

When extended sources other than the cluster were detected in the images, we attempted to determine if they were physically associated with the cluster or were unrelated projected sources. We did this by querying the NASA Extragalactic Database⁹ to see if the source has a known counterpart.

• A773 – For ObsID 5006, due to flaring we only use the first 15 ks of the observation.

• MS0906 – At 9:08:58.02, 11:01:58.36 (RA, Dec = 137.24, 11.03), there is a large region of diffuse extended emission (which corresponds to cluster A750 Abell et al. 1989) which we mask with a 230″ radius circle.

• ZW2701 – At 9:53:05.85, 51:49:16.39 (148.27, 51.82) there is a region of extended emission that we mask with a 150″ radius circle.

• A963 – For ObsID 903, a temperature of 0.28 keV was used when fitting the APEC model to the soft background residuals. We note that the standard temperature used for this is 0.18 keV, as it is found that this temperature generally models these residuals well (Giles et al. 2015). The higher temperature implies that this observation includes a hotter region of the galactic foreground emission.

• ZW3146 – For ObsID 9371, due to flaring we only use the first 34 ks of the observation.

• A1423 – For ObsID 538, due to flaring we only use the first 11 ks of the observation.

• A1553 – For ObsID 12254, due to flaring we only use the first 10 ks of the observation.

• A1682 – At 13:06:59.904, 46:31:40.65 (196.75, 46.53) there is a cluster galaxy (Morrison et al. 2003) and at 13:07:13.47, 46:29:02.31 (196.81, 46.48) there is another small region of extended emission. Both sources are masked by radius 40″ circles.

• A1763 – At 13:34:52.8, 40:57:21.6 (203.72, 40.96) we use a 90″ radius circle to mask a region of extended emission, which is likely to be associated with a previously known X-ray source (Evans et al. 2010).

• A1930 – At 14:32:42.72, 31:33:50.4 (218.18, 31.56) there is a region of extended emission that we mask with a circle of radius 160″. The extended emission is the cluster RM J143242.6+313407.1 (Rozo et al. 2015) at redshift z = 0.137, and is likely part of the same dark matter halo as A1930 (z = 0.1308).

• A2009 – For ObsID 10438, a temperature of 0.27 keV was used when fitting the APEC model to the soft background residuals.

• A2069 – There is a separate cluster at 15:24:25.846, +30:00:16.039 (231.11, 30.00) at redshift z = 0.119 (MaxBCG J231.10029+30.00604, Koester et al. 2007); it is likely part of the same dark matter halo as A2069. We mask it with a rectangle of 23′ by 14′ at an inclination angle of 35°. There is either a group of point sources or a filament at 15:23:38.87, +29:58:30.226 (230.91, 29.98) that we mask using an ellipse. In addition, for ObsID 4965, due to flaring we only use the first 40 ks of the observation.

• A2261 – At 17:22:12.78, 32:06:36.95 (260.55, 32.11) there is a region of diffuse extended emission that we mask with a 80″ radius circle. This source is associated with a galaxy cluster (GMBCG J260.55436+32.11438) at z = 0.304 (Hao et al. 2010), so is not associated with A2261 which is at redshift z = 0.2242.

• RXJ1720 – For ObsID 4361, due to flaring we only use the first 15 ks of the observation.

Appendix B: Temperature and density profiles

In this Appendix we expand on the method used to model the temperature and density profiles used to determine the hydrostatic masses in Sect. 3.

B.1. Temperature Profiles

The model used for the 3D distribution of the temperature (Vikhlinin et al. 2006) is sufficient to describe a temperature decline (if present) in the central core region of the cluster in addition to the rest of the cluster temperature profile outside the central cooling region of a cluster. The part outside the central cooling region is modelled by

$$\begin{array}{c}\hfill t(r)=-\frac{{(r/r_t)}^{-a}}{{[1+{(r/r_t)}^b]}^{c/b}},\end{array}$$(B.1)

where r is the radius of the cluster and all other symbols are free parameters: a models the slope at small radii, c models the slope at large radii, and b models the width of the transition between these slopes. This transition occurs at a radius r_t. The decrease in temperature found in the central cooling region of many clusters can be modelled as

$$\begin{array}{c}\hfill t_{\mathit{cool}}(r)=-\frac{(x+T_{\mathit{min}}/T_0)}{x+1},x={\left(\frac{r}{r_{\mathit{cool}}}\right)}^{a_{\mathit{cool}}},\end{array}$$(B.2)

where r is cluster radius and all other symbols are free parameters: T₀ is the normalisation, T_min is the minimum temperature, r_cool models radius of the central cooling region and a_cool models the slope of the cool region which extends to a radius r_cool.

The 3D temperature profile model is a product of equations B.1 and B.2

$$\begin{array}{c}\hfill T_{3D}(r)=T_0{t}_{\mathit{cool}}(r)t(r).\end{array}$$(B.3)

B.2. ICM density profiles

The model used to describe the 3D distribution of the density (Vikhlinin et al. 2006) is given by:

$$\begin{array}{c}\hfill n_{\mathrm{e}}{n}_{\mathrm{p}}=n_0^2\frac{{(r/r_c)}^{-\alpha }}{{(1+r^2/r_c^2)}^{3\beta -\alpha /2}}\frac{1}{{(1+r^{\gamma }/r_s)}^{ϵ/\gamma }}+\frac{n_{02}^2}{{(1+r^2/r_{c2}^2)}^{3{\beta }_2}},\end{array}$$(B.4)

with r being the cluster radius and the nine free parameters able to model both the central and outer parts of the density profile sufficiently. We fix γ = 3 as in Vikhlinin et al. (2006), and all the other parameters are free. n_e and n_p are the electron and proton density respectively, and r is the distance from the centre of the cluster. The first term in the equation is a modification of the β-model (Cavaliere & Fusco-Femiano 1978) that allows for a power-law type cusp, with slope α at small radii, which is expected in the centres of clusters that are dynamically relaxed. The second term models the change of the slope by ϵ near the radius r_s, and the width of this transition region is controlled by γ. The final term is another β-model component, which increases the freedom of the overall model near the cluster centre.

B.3. Projection of 3D profiles

The 3D temperature and density profiles must be projected before they can be compared with the observed profiles. For this projection we assume spherical symmetry (see e.g. Ettori et al. 2013, for the geometry of this projection). We note that whenever we use an abundance profile in our analysis, we use the measured projected abundance profile (we lack the very high quality data needed to deproject abundance profiles), and at large radii set the abundance to 0.3 Z_⊙, in line with observed flattening of abundance profiles at large radii (Leccardi & Molendi 2008).

The projected temperature profile was computed using the method described in Vikhlinin et al. (2006). Briefly, the temperature that would be measured when fitting a single temperature thermal plasma model to the projected emission is computed by weighting the different temperature components of the 3D model along the line of sight by their emissivity. The emissivity is derived from the projected abundance and 3D density model, and the weighting algorithm of Vikhlinin (2006) was used.

The 3D density profile is projected to produce a model surface brightness profile. The count rate from each volume element along the line of sight is computed assuming an APEC plasma emission model with the normalisation set by the density in the volume element, and the temperature given by the 3D temperature profile model. This conversion from 3D density to projected surface brightness uses the appropriate instrument response for the detector location.

In the fitting process, the projection of the temperature and density profiles are performed simultaneously for a given set of model parameters to produce projected temperature and surface brightness profiles that are compared with the data. This is an improvement over the method used in Maughan et al. (2016), where a normalised average temperature profile was assumed when converting from density to emissivity.

Appendix C: Mass profile plots

Figure C.1 shows the X-ray hydrostatic and caustic mass profiles for the 44 clusters in our sample.

Fig. C.1. The caustic and hydrostatic mass profiles for each cluster as black solid and green dashed lines respectively in the upper panels of each plot. The ratio of the hydrostatic to caustic mass, MX/MC, are shown in the lower panels of each plot. 1σ uncertainties are shown by the shaded regions. The vertical black line is at the value of R500 as calculated from the hydrostatic mass profile; the solid red vertical line is at the outer radius of the measured temperature profile (note that hydrostatic masses beyond this radius are based on extrapolation).

Fig. C.1. continued.

Fig. C.1. continued.

Fig. C.1. continued.

Fig. C.1. continued.

Fig. C.1. continued.

All Tables

All Figures

Fig. 1. Graphical model summarising the mass modelling framework that we present in Sect. 3.3. The symbols are also described in that section; in summary, the symbol μ represents mass, κ is the bias term, δ is the scatter term and σ is the measurement error term (all in log space), and the C and X subscripts denote whether the quantity refers to that from the hydrostatic or caustic method. The refers to an observed quantity. 𝒩 refers to a normal distribution, and ‘∼’ denotes ‘is distributed as’.

In the text

	Fig. 2. The hydrostatic versus the caustic masses are plotted for the full sample of clusters (with masses calculated at the hydrostatic R500). NRCC clusters are red circles and RCC clusters are blue triangles. The solid black line shows the best fitting model to the full data, with the shaded region indicating the 1σ uncertainty. The 1:1 line is also plotted as a dashed black line.
In the text

	Fig. 3. The MX/MC profiles for all 44 clusters in our sample, scaled in radius by the hydrostatic R500. NRCC clusters are solid red lines and RCC clusters are blue dot dashed lines. The 1:1 line is also plotted as a thin dashed black line. The thick dashed black line shows the average MX/MC ratio, and the shaded region shows the 1σ uncertainty.
In the text

Fig. 4. We show the M500, X/M500, C ratio for all 44 clusters in our sample. NRCC clusters are red circles and RCC clusters are blue triangles. The solid black line and shaded region show the average M500, X/M500, C ratio at the hydrostatic R500. 1σ errors are shown. The clusters are ranked by Ngal with the highest at the top, and the horizontal lines separate the clusters into the low, mid and high Ngal subsets (see Sect. 4.2).

In the text

	Fig. 5. The hydrostatic and caustic masses measured within the hydrostatic R500 are plotted for the full sample. This is the same data as plotted in Fig. 2, but now showing the Ngal subsets. The high, mid and low Ngal bins are green squares, blue triangles and red circles respectively. The 1:1 line is also plotted as a dashed black line.
In the text

	Fig. 6. Left: fgas measured at 1 Mpc using caustic mass in the denominator is plotted against Ngal. Right: same plot, but with fgas calculated with hydrostatic mass in the denominator. In both panels, the black dashed lines indicate fgas = 0.10 and fgas = 0.15 for reference, encompassing the approximate expected range of fgas.
In the text

Fig. 7. The ratio of the caustic masses (measured at the hydrostatic R500) from Sohn et al. (2020) to those from Rines et al. (2013) are plotted as a function of the difference in Ngal between the analyses. The clusters are split into subsets containing those where the minimum Ngal used in the two analyses was at least 200 (green squares) and those where the minimum Ngal was less then 120 (red circles). The lines show the best fits to each subset.

In the text

	Fig. 8. The average MX/MC profiles for the high Ngal, mid Ngal and low Ngal clusters in green lines, blue lines, and red dot-dashed lines respectively. The 1:1 line is also plotted as a thin black dashed line. The thick dashed black line shows the average MX/MC profile for all clusters. The profiles are scaled in radius by the hydrostatic R500, and the shaded regions show the 1σ uncertainties.
In the text

Fig. 9. We show the hydrostatic masses versus the caustic masses for the 14 clusters in our sample with the highest Ngal values, both calculated at the hydrostatic R500. We show RCC clusters as blue triangles and NRCC clusters as red circles. The solid black line shows the best fitting model, with the shaded region indicating the 1σ uncertainty. The 1:1 line is also plotted as a dashed black line.

In the text

Fig. C.1. The caustic and hydrostatic mass profiles for each cluster as black solid and green dashed lines respectively in the upper panels of each plot. The ratio of the hydrostatic to caustic mass, MX/MC, are shown in the lower panels of each plot. 1σ uncertainties are shown by the shaded regions. The vertical black line is at the value of R500 as calculated from the hydrostatic mass profile; the solid red vertical line is at the outer radius of the measured temperature profile (note that hydrostatic masses beyond this radius are based on extrapolation).

In the text

	Fig. C.1. continued.
In the text

	Fig. C.1. continued.
In the text

	Fig. C.1. continued.
In the text

	Fig. C.1. continued.
In the text

	Fig. C.1. continued.
In the text