Debiased Galaxy Cluster Pressure Profiles from X-Ray Observations and Simulations

For a spherically symmetrical cluster, hydrostatic equilibrium occurs when the the force of gravity exerted on gas in the cluster is balanced by the gradient of gas pressure:

$$\begin{eqnarray}&&\displaystyle \frac{{dP}(r)}{{dr}}=-{\rho }_{\mathrm{gas}}(r)\displaystyle \frac{{GM}(\lt r)}{{r}^{2}},\end{eqnarray} \tag{ 1 }$$

with the gravitational constant, G, enclosed mass profile, M(<r), gas pressure profile, P(r), and gas density profile, ρ_gas(r), all with respect to the distance r from the center of the cluster. A combination of X-ray observations like XMM-Newton, CHANDRA, and analysis techniques taking into account projection and point-spread function effects has achieved high-resolution measurements of the radial electron density profiles, n_e(r), and the radial temperature profiles, T(r), of galaxy clusters (e.g., Vikhlinin et al. 2006; Croston et al. 2008), which can be used to determine the radial electron pressure profile, P_e(r), by assuming an ideal gas equation of state, P_e(r) = k_B n_e(r)T(r). Given the electron pressure, the gas thermal pressure P_th is defined by P_th(r) = P_e(r)μ_e/μ where μ is the mean mass per gas particle, and μ_e is the mean mass per electron.

In addition to the thermal motion of the gas, other sources of gas pressure—including viralized bulk motion, turbulence, cosmic rays, and magnetic fields—also provide non-trivial pressure support (e.g., Ensslin et al. 1997; Churazov et al. 2008; Brüggen & Vazza 2015). For realistic equilibrium systems, the gas pressure, P, in Equation (1) is replaced by P = P_th + P_nth, where P_nth refers to any non-thermal pressure acting on the intracluster gas. X-ray-measured cluster masses are derived from the assumption of hydrostatic equilibrium with only thermal gas pressure, which means that the contribution of non-thermal pressure can cause X-ray measurements to underestimate cluster masses systematically.

Numerical simulations provide a vital resource for characterizing the mass bias as the properties of simulated galaxy clusters are known exactly. Barnes et al. (2020) developed the Mock-X analysis framework, which can generate synthetic X-ray images and derives halo properties (e.g., gas density and temperature profiles) via observational methods, which can be used to derive hydrostatic mass in mock X-ray observations. Hydrostatic mass bias is equal to the ratio of the hydrostatic mass to the "true" (overdensity) mass of simulated clusters identified through SUBFIND (e.g., Springel et al. 2001; Dolag et al. 2009) in simulations. Studies (e.g., Henson et al. 2017; Barnes et al. 2020) also point out that the bias of hydrostatic mass estimated with density and temperature profiles derived from the spectroscopic analysis show a much stronger mass-dependence than those estimated from the true mass-weighted temperature profiles. These simulated spectroscopic temperatures emulate the observational procedure for measuring X-ray temperatures and thus we compare against them in this analysis.

A number of the aforementioned numerical studies have measured ${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$/${M}_{500{\rm{c}}}^{\mathrm{True}}.$ However, observationally, one only has access to ${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}.$ This means that we must invert these relations to give ${M}_{500{\rm{c}}}^{\mathrm{True}}/{M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$ as a function of ${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$ (a deceptively complex task). Here, M_500c is the "overdensity mass" and corresponds to the mass within a spherical boundary which has an average density equal to 500 times the critical density, ρ_crit.

In this work, we present an efficient approach to estimate the true cluster mass by utilizing both the X-ray and "true" masses of simulated clusters, ${M}_{500{\rm{c}}}^{\mathrm{True}}$ from Barnes et al. (2020). We adopt a power-law model for the scaling relation between ${M}_{500{\rm{c}}}^{\mathrm{True}}$ and ${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$. This is a linear model in logarithmic scale. For convenience, we denote

$$\begin{eqnarray}\begin{array}{rcl}{S}_{{\rm{X}}} & = & \mathrm{ln}({M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}/{M}_{0}),\\ {S}_{{\rm{M}}} & = & \mathrm{ln}({M}_{500{\rm{c}}}^{\mathrm{True}}/{M}_{0}),\end{array}\end{eqnarray} \tag{ 2 }$$

with M₀ = 3 × 10¹⁴ M_⊙. For fixed S_X, we assume a linear relation between S_M and S_X where the error, , follows a Gaussian distribution:

$$\begin{eqnarray}&&{S}_{{\rm{M}}}={{aS}}_{{\rm{X}}}+b+\epsilon ,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\epsilon \sim \mathrm{Norm}(0,{\sigma }^{2}),\end{eqnarray} \tag{ 4 }$$

where a, b, σ are free parameters.

We notice that clusters drawn from simulations are selected in terms of a certain mass threshold, which means we also need to consider this selection effect in our model when fitted to simulation data. For given S_X and S_M of a simulated cluster, we use a truncated normal distribution to model the likelihood

$$\begin{eqnarray}&&p({S}_{{\rm{M}}}| {S}_{{\rm{X}}},{S}_{{\rm{T}}},{\boldsymbol{\theta }})=\displaystyle \frac{A({S}_{{\rm{X}}},{S}_{{\rm{T}}},\theta )}{\sqrt{2\pi {\sigma }^{2}}}\exp \left[\displaystyle \frac{{\left({{aS}}_{{\rm{X}}}+b-{S}_{{\rm{M}}}\right)}^{2}}{2{\sigma }^{2}}\right],\end{eqnarray} \tag{ 5 }$$

where θ = (a, b, σ) denotes the free parameters. S_T is the truncation parameter defined by ${S}_{{\rm{T}}}=\mathrm{log}({M}_{500{\rm{c}}}^{{\rm{T}}}/{M}_{0})$ and ${M}_{500{\rm{c}}}^{{\rm{T}}}$ is the mass threshold for a given simulated cluster sample. A(S_T, S_X, θ ) is the normalization factor for a normal distribution, Norm(aS_X + b, σ²), truncated with a lower bound S_T:

$$\begin{eqnarray}&&A({S}_{{\rm{T}}},{S}_{{\rm{M}}},{\boldsymbol{\theta }})={\left[1-{\rm{\Phi }}\left(\displaystyle \frac{{S}_{{\rm{T}}}-{{aS}}_{{\rm{X}}}-b}{\sigma }\right)\right]}^{-1},\end{eqnarray} \tag{ 6 }$$

where Φ denotes the cumulative distribution function of standard normal. We set

$$\begin{eqnarray}&&a,b\sim { \mathcal U }(-5,5)\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\sigma \sim { \mathcal U }(0,5)\end{eqnarray} \tag{ 8 }$$

as priors, where ${ \mathcal U }$ denotes the uniform distribution. We can then write out the posterior for the parameters:

$$\begin{eqnarray}&&p({\boldsymbol{\theta }}| D)\propto p(a)p(b)p(\sigma )\displaystyle \prod _{\alpha }\left[\displaystyle \prod _{i=1}^{{N}_{\alpha }}p({S}_{{\rm{M}},\alpha }^{i}| {S}_{{\rm{X}},\alpha }^{i},{S}_{{\rm{T}},\alpha },{\boldsymbol{\theta }})\right],\end{eqnarray} \tag{ 9 }$$

where D = {D₁, D₂, D₃} is a data vector of log-scaled masses of simulated cluster samples drawn from IllustrisTNG, BAHAMAS, and MACSIS simulations denoted by α = 1, 2, 3, and ${D}_{\alpha }=\{({S}_{{\rm{M}},\alpha }^{i},{S}_{{\rm{X}},\alpha }^{i}),i=1,...\,{N}_{\alpha }\}$. Each simulation uses a different S_T,α . Mass thresholds, ${M}_{500{\rm{c}}}^{{\rm{T}}},$ are set to be 10¹⁴ M_⊙ for IllustrisTNG and BAHAMAS, and 4 × 10¹⁴ M_⊙ for MACSIS. Details about these simulations can be found in Barnes et al. (2020).

We note that the IllustrisTNG, BAHAMAS, and MACSIS simulations adopt different numerical methods or subgrid physics, which may introduce differences in the derived cluster profiles and systematics in the mass estimation of the mock X-ray observation. This will be accounted for in the intrinsic scatter in our regression model for the relation of ${M}_{500{\rm{c}}}^{\mathrm{True}}$ versus ${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$ since the fit is performed on all simulations simultaneously.

We explore the parameter space by Markov chain Monte Carlo, using emcee (Foreman-Mackey et al. 2013) for the sampling. We discard the initial steps suggested by the integrated autocorrelation time (Foreman-Mackey et al. 2019), which estimates the number of steps that are needed before the chain "forgets" where it started. This step ensures the samples are well "burnt-in." Regression results for the linear model and uncertainty are reported in Table 1.

A small fraction, ∼9%, of simulated clusters have abnormally high or abnormally low ${M}_{500{\rm{c}}}^{\mathrm{True}}$/${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$, suggesting ratios outside the observed range (e.g., Miyatake et al. 2019). Most cases appear in low-mass clusters, and may be due to the numerical noise when resolving the X-ray mass of simulated clusters from synthetic images. The steep slope of the mass function causes these unreliable low-mass data points to significantly influence the mean ${M}_{500{\rm{c}}}^{\mathrm{True}}$ at high ${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$. In addition to the numerical noise, merging events, or certain active galactic nucleus activities in the unrelaxed clusters could lead to a less spherical cluster. The thermodynamic profiles and corresponding X-ray masses of these less regular clusters will be recovered with more systematic uncertainty because clusters' profiles and masses are derived assuming spherical symmetry in the Mock-X analysis (Barnes et al. 2020), which could also result in an extreme value of ${M}_{500{\rm{c}}}^{\mathrm{True}}$/${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$. For a more concrete conclusion, detailed studies are required for these peculiar cluster samples with abnormal values of ${M}_{500{\rm{c}}}^{\mathrm{True}}$/${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$ in future work.

To mitigate the effect introduced by the simulated clusters with extreme values of ${M}_{500{\rm{c}}}^{\mathrm{True}}$/${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$, we iteratively remove outlier clusters falling outside the 2σ region of the regression results until the prediction for ${M}_{500c}^{\mathrm{True}}$ derived from the linear model for ${M}_{500c}^{\mathrm{True}}$ versus ${M}_{500c}^{{\rm{X}} \mbox{-} \mathrm{ray}}$ converges to 1% agreement with the previous iteration. This is performed for clusters within the X-ray mass range ${M}_{500c}^{{\rm{X}} \mbox{-} \mathrm{ray}}={10}^{14}\mbox{--}{10}^{15}{M}_{\odot }$. To test the impact of this method for removing outliers, we also use a truncated t-distribution (e.g., Pfanzagl & Sheynin 1996) to model the uncertainty, , which is another approach to alleviate the effect of outlier samples.

In Figure 1, we plot S_M versus S_X for the IllustrisTNG, BAHAMAS, and MACSIS cluster samples. We also show the regression results for the linear relation between log-scaled "true" and X-ray masses, considering both truncation effects and the influence of outlier clusters. The intrinsic scatter in our linear model is determined by the parameter σ. We find the slope parameter a = 1.079 is greater than 1, which indicates the ratio of ${M}_{500{\rm{c}}}^{\mathrm{True}}$ to ${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$ is mass-dependent, and hydrostatic mass bias increases with cluster mass. We also plot the regression results for the linear relation by modeling the uncertainty, , with a truncated t-distribution. For comparison, we also show regression results without removing outlier clusters and without modeling truncation effects.

Zoom In Zoom Out Reset image size

The alternative fit for ${M}_{500{\rm{c}}}^{\mathrm{True}}$ versus ${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$, which does not perform clipping but uses a truncated t-distribution to account for outliers, is in good agreement with the results of iterative 2σ clipping methods. Regression results of the two methods find similar values for the slope parameter and the discrepancy between ${M}_{500{\rm{c}}}^{\mathrm{True}}$ for ${10}^{14}{M}_{\odot }\lt {M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}\lt {10}^{15}{M}_{\odot }$ is at the ∼ 1% level. Comparing with the fit which did not account for outliers, we find outlier removal has a modest but statistically significant effect on fit results. We also find that failing to account for mass selection effects results in a bad fit to the simulation data and a substantially different power-law index.

When we fit an analytical model like a GNFW profile to the radial pressure profile of a galaxy cluster, a common approach taken is to normalize the pressure and radius by the characteristic pressure P_500c and radius R_500c , both of which can be directly computed at a given cluster mass. If the mass of galaxy clusters in X-ray measurements suffers from hydrostatic bias, the characteristic pressure and radius will as well. For convenience, we define a new variable for the hydrostatic mass bias,

$$\begin{eqnarray}&&{B}_{{\rm{M}}}={M}_{500{\rm{c}}}^{\mathrm{True}}/{M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}},\end{eqnarray} \tag{ 10 }$$

then the radius bias, B_R, and pressure bias, B_P, can be obtained from scaling relations, although the latter relies on the assumption of a specific model for the pressure profile. For a spherical cluster, ${R}_{500{\rm{c}}}\propto {M}_{500{\rm{c}}}^{1/3}$, so hydrostatic bias for cluster radius is defined by

$$\begin{eqnarray}&&{B}_{{\rm{R}}}={B}_{{\rm{M}}}^{1/3}.\end{eqnarray} \tag{ 11 }$$

If we assume that pressure follows a GNFW profile given by $P(r)={P}_{500{\rm{c}}}({M}_{500{\rm{c}}},z){\mathbb{P}}(x)$ (Nagai et al. 2007a), where x = r/R_500c and

$$\begin{eqnarray}\begin{array}{rcl}{P}_{500{\rm{c}}}({M}_{500{\rm{c}}},z) & = & 1.65\times {10}^{-3}h{\left(z\right)}^{8/3}\\ & & \times \,{\left[\displaystyle \frac{{M}_{500c}}{3\times {10}^{14}{M}_{\odot }}\right]}^{2/3}{h}_{70}^{2}\,\mathrm{keV}\,{\mathrm{cm}}^{-3},\end{array}\end{eqnarray} \tag{ 12 }$$

${\mathbb{P}}(x)$ is the scaled profile, with the form

$$\begin{eqnarray}&&{\mathbb{P}}(x)=\displaystyle \frac{{P}_{0}}{{\left({c}_{500}x\right)}^{\gamma }{\left[1+{\left({c}_{500}x\right)}^{\alpha }\right]}^{\left(\beta -\gamma \right)/\alpha }},\end{eqnarray} \tag{ 13 }$$

where c₅₀₀ is the concentration, P₀ is the normalization parameter, and the parameters α, β, γ determine the power-law slopes of different region of the cluster. Since ${P}_{500{\rm{c}}}\propto {M}_{500{\rm{c}}}^{2/3}$ according to Equation (12), the pressure bias is

$$\begin{eqnarray}&&{B}_{{\rm{P}}}={B}_{{\rm{M}}}^{2/3}.\end{eqnarray} \tag{ 14 }$$

The bias parameters B_R and B_P can be used to debias GNFW fits to X-ray measurements of thermal pressure profiles by rescaling c₅₀₀ and P₀ with the following bias correction factors:

$$\begin{eqnarray}&&{c}_{500}={c}_{500}^{\mathrm{bias}}\times {B}_{{\rm{R}}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{P}_{0}={P}_{0}^{\mathrm{bias}}/{B}_{{\rm{P}}}.\end{eqnarray} \tag{ 16 }$$

We note that radius and pressure biases have a one-to-one relation with B_M, so the uncertainty in B_M can be converted to B_R and B_P by

$$\begin{eqnarray}&&{\sigma }_{\mathrm{ln}{B}_{{\rm{R}}}}={\sigma }_{\mathrm{ln}{B}_{{\rm{M}}}}/3,{\sigma }_{\mathrm{ln}{B}_{{\rm{P}}}}=2{\sigma }_{\mathrm{ln}{B}_{{\rm{M}}}}/3.\end{eqnarray} \tag{ 17 }$$

We apply our linear model for S_M versus S_X to the hydrostatic X-ray masses of the REXCESS cluster sample to estimate the true masses of these clusters. REXCESS is a representative sample of local clusters at redshifts 0.0 < z < 0.2 which spans a mass range of 10¹⁴ M_⊙ < M_500c < 10¹⁵ M_⊙ (Arnaud et al. 2010). REXCESS clusters are drawn from the REFLEX catalog and were studied in-depth by the XMM-Newton Large Programme. A description of the REFLEX sample and of XMM-Newton observation details can be found in Böhringer et al. (2007). We correct the hydrostatic mass of 31 local clusters from the REXCESS sample measured by X-ray observation:

$$\begin{eqnarray}&&{M}_{500{\rm{c}}}^{\mathrm{True}}/{M}_{0}={e}^{b}\times {\left({M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}/{M}_{0}\right)}^{a},\end{eqnarray} \tag{ 18 }$$

where a and b are the regression parameters for the linear model reported in Table 1. We find that the X-ray-measured hydrostatic masses of clusters in the REXCESS sample are underestimated by approximately 7% on average. The bias climbs from 0% to 15% as cluster X-ray mass increases from 10¹⁴ M_⊙ to 10¹⁵ M_⊙.

The regression parameter σ is the intrinsic scatter in ${S}_{{\rm{M}}}=\mathrm{ln}({M}_{500{\rm{c}}}^{\mathrm{True}}/{M}_{0})$ and can be used to characterize the uncertainty in the corrected mass of REXCESS clusters at a given predicted ${M}_{500{\rm{c}}}^{\mathrm{True}}:$

$$\begin{eqnarray}&&{\sigma }_{{M}_{500{\rm{c}}}^{\mathrm{True}}}\simeq {M}_{500{\rm{c}}}^{\mathrm{True}}\sigma .\end{eqnarray} \tag{ 19 }$$

With the first-order approximation, this scatter yields significant uncertainties for individual objects, around ≈20%, for corrected cluster masses. σ also defines scatter in the mass bias B_M, radius bias B_R, and pressure bias B_P in log scale:

$$\begin{eqnarray}&&{\sigma }_{\mathrm{ln}{B}_{{\rm{M}}}}=\sigma ,{\sigma }_{\mathrm{ln}{B}_{{\rm{R}}}}=\sigma /3,{\sigma }_{\mathrm{ln}{B}_{{\rm{P}}}}=2\sigma /3.\end{eqnarray} \tag{ 20 }$$

These allow us to estimate the modeling uncertainty in the debiased pressure, radius, and mass.

The UPP is a model for ICM thermal pressure profiles developed by Arnaud et al. (2010) which was calibrated off the REXCESS sample. For each cluster in the sample, the pressure profile—derived along with the X-ray measurements of gas density and temperature profiles—is scaled with the characteristic pressure P_500c and cluster radius R_500c. As discussed in Section 2.2, both R_500c and P_500c are dimensional rescalings of M_500c, which itself is measured from a M_500c − Y_X relation (Kravtsov et al. 2006; Arnaud et al. 2007; Nagai et al. 2007a) which was calibrated on biased hydrostatic mass estimates. Note that Arnaud et al. (2007) do not use the REXCESS sample, which could potentially allow selection bias to creep in. The UPP model is widely used for characterizing cluster masses in SZ surveys (e.g., Hasselfield et al. 2013; Planck Collaboration et al. 2016a; Hilton et al. 2018) and is expressed as

$$\begin{eqnarray}\begin{array}{rcl}P(x,{M}_{500{\rm{c}}},z) & = & {P}_{500{\rm{c}}}({M}_{500{\rm{c}}},z)\\ & & \times \,{\mathbb{P}}(x){\left[\displaystyle \frac{{M}_{500{\rm{c}}}}{3\times {10}^{14}{h}_{70}^{-1}{M}_{\odot }}\right]}^{{\alpha }_{{\rm{P}}}(x)},\end{array}\end{eqnarray} \tag{ 21 }$$

with variables taking the same meaning as in Section 2.2. The empirical term, ${({M}_{500{\rm{c}}}/3\times {10}^{14}{h}_{70}^{-1}{M}_{\odot })}^{{\alpha }_{{\rm{P}}}(x)}$, reflects the deviation from standard self-similar scaling with α_P(x) = 0.22/(1 + 8x³). A GNFW profile, ${\mathbb{P}}(x)$, is fit against the (geometric) mean profile of the scaled REXCESS sample.

The hydrostatic bias that we found for M_500c in the REXCESS sample is transferred to the normalization of observed pressure profiles through the resultant changes in P_500c and R_500c. For each REXCESS cluster, we use the GNFW pressure profile provided in Arnaud et al. (2010) and rescale P₀, and c₅₀₀ according to Equations (15) and (16) to get the debiased fits for each cluster. We then evaluate the geometric mean of the scaled profiles, P_m, and fit it with a GNFW model in the log–log plane. We also estimate the uncertainty in the mean profile by approximating the uncertainty in each corrected pressure profile via lognormal distributions with variances ${\sigma }_{\mathrm{ln}R}$ and ${\sigma }_{\mathrm{ln}P}$. Moreover, we use this uncertainty to define the 68% range for the mean profile confined by a high profile, P_h, and a low profile, P_l.

We fit new GNFW models to the mean, high, and low profiles discussed above, fixing the outer slope parameter to β = 5.490 as was done in the original UPP model. In Arnaud et al. (2010), the GNFW model of the UPP is fitted to the observed average scaled profile in the radial range [0.031]R_500c, combined with the average simulation profile beyond R_500c which is crucial for determining the outer slope β. In our paper, the GNFW model is fitted to the debiased observed profiles within R_500c, but we lack information beyond this radius. So we choose to keep β as same as its original value in the UPP model. The best-fitting parameters of the GNFW models for P_m, P_h, and P_l are reported in Table 2.

In the top panel of Figure 2, we plot corrected GNFW fits to the DPPs for each of the 31 REXCESS clusters. As discussed in Section 3.1, the scatter in B_M is significant and introduces non-negligible uncertainty to the DPPs of the REXCESS sample. We also show the uncertainty in the DPP for each RECXESS cluster considering the uncertainty as determined by ${\sigma }_{\mathrm{ln}{B}_{{\rm{R}}}}$ and ${\sigma }_{\mathrm{ln}{B}_{{\rm{P}}}}$. We show the geometric mean of these scaled profiles, the fit to this curve, and the UPP model for comparison. The dispersion in these scaled pressure profiles is significant in the core both before and after debiasing regions due to the various dynamical states, including both the cool core and morphologically disturbed clusters of the REXCESS sample (Arnaud et al. 2010). The mean of the debiased scaled pressure profiles and its GNFW fit, ${\mathbb{P}}(x)$, is lower than in the original UPP model. In the bottom panel of Figure 2, we plot the fractional difference between both the UPP and the debiased mean scaled profile against our best fit to the latter. We also show uncertainty in the pressure model with the red semitransparent region.

Zoom In Zoom Out Reset image size

The UPP is ≈5% higher than the mean of the DPP in the center of the cluster and gradually climbs to 20% at R_500c, and reaches almost ≈35% at the outermost outskirts. Only weak scattering is found for the adjusted scaled pressure profile compared to the uncertainty of scaled pressure profile of each REXCESS cluster, which is due to the assumption of using Gaussian approximating the uncertainty of individual profile in logarithmic scale at fixed radii, and uncertainty of the mean decreases with the growth of the sample size of the REXCESS clusters.

We also explore whether the REXCESS pressure profiles deviate from self-similarity by studying their radial variation as a function of mass. To do this, we look for mass trends in P(x)/P_500c in our debiased profiles. We evaluate these profiles at x = r/R_500c = 0.1, 0.2, 0.4, 0.8. This range of radii avoids either too small or too large values of x. We avoid larger scaled radii because X-ray measurements pressure profiles in REXCESS clusters rarely get beyond R_500c (Arnaud et al. 2010). ⁵ We avoid taking a smaller value of x because the REXCESS sample contains systems with various dynamical states which can alter the state of gas in the center of the cluster. Following Arnaud et al. (2010), we fit a power-law of the form $P/{P}_{500{\rm{c}}}\propto {M}_{500{\rm{c}}}^{{\alpha }_{{\rm{P}}}(x)}$ to each set of points weighted by uncertainties on both cluster masses and pressure following the orthogonal regression approach, proposed for the analysis when both the dependent and the independent variables are random. The best-fitting results and comparison with the UPP model are represented in Table 3.

In Figure 3, we show the results of this fit, with different colors representing different scaled radii and error bars representing uncertainty due to the intrinsic scatter in B_M. We show the best-fit power-laws to each set of points and the values of their power indices.

Zoom In Zoom Out Reset image size

Our study of the debiased scaled pressure profiles of the REXCESS cluster sample finds that α_P(x) at all radii are consistent with zero, which means a less significant deviation from standard self-similarity compared to the UPP model. We can observe a radial dependence of α_P(x) similar to that found in the UPP model. However, this term in UPP is treated as a second-order deviation term in addition to a constant modification of the standard self-similarity, α_P ∼ 0.12, which can be neglected in the first-order approximation. Based on the discussion above, we see no evidence for deviations from self-similarity, which would require the mass-dependent term in Equation (21). We modify the UPP by eliminating the deviation term and get a simplified model for ICM pressure profiles, the DPP:

$$\begin{eqnarray}&&{P}_{\mathrm{DPP}}(x,{M}_{500{\rm{c}}},z)={P}_{500{\rm{c}}}({M}_{500{\rm{c}}},z)\times {\mathbb{P}}(x).\end{eqnarray} \tag{ 22 }$$

Here, parameters take on the same meaning as in P_m in Table 2.

The spherical volume-integrated Compton parameter, Y_sph, of a cluster is the integral of the gas's thermal pressure profile over a spherical region and is defined as

$$\begin{eqnarray}&&{Y}_{\mathrm{sph}}(R)=\displaystyle \frac{{\sigma }_{{\rm{T}}}}{{m}_{{\rm{e}}}{c}^{2}}{\int }_{0}^{R}4\pi {P}_{{\rm{e}}}(r){r}^{2}{dr},\end{eqnarray} \tag{ 23 }$$

where σ_T is the Thomson cross-section, m_e is the electron mass, and P_e is the thermal electron pressure. Since the pressure is directly related to the depth of cluster gravitation potential, Y_sph is closely related to the mass of the cluster. Studies (e.g., Da Silva et al. 2004; Nagai 2006) find a low intrinsic scatter in the relation between integrated Compton parameter and cluster mass, indicating that Y_sph serves as a good proxy for cluster mass. The Y_sph–M relation was previously modeled with a power law (Kravtsov et al. 2006; Arnaud et al. 2007; Nagai et al. 2007a). Accordingly, we parameterize the Y_sph(R_500c) − M_500c relation as

$$\begin{eqnarray}&&h{\left(z\right)}^{-2/3}{Y}_{\mathrm{sph}}({R}_{500{\rm{c}}})={10}^{A}{\left[\displaystyle \frac{{M}_{500{\rm{c}}}}{3\times {10}^{14}{M}_{\odot }}\right]}^{\alpha }{h}_{70}^{-5/2}{\mathrm{Mpc}}^{2}.\end{eqnarray} \tag{ 24 }$$

We fit Equation (24) to the X-ray-measured Compton parameter and the biased X-ray hydrostatic masses of the REXCESS sample and find that α = 1.790 ± 0.015 and A = −4.739 ± 0.003. The Y_sph–M relation can be derived from the UPP model by combining Equation (21) for the UPP and Equation (23) and gives α = 1.787, and A = −4.745. The analytical calculations based on the UPP and direct fits to observation data are in excellent agreement: both claim a deviation from the slope predicted by self-similarity, α_s = 5/3, of approximately Δα = α − α_s ≈ 0.12. Notice this deviation Δα corresponds to the α_P (x) for the pressure model, which is characterized by a function of cluster mass and radius; however, Arnaud et al. (2010) showed this term can be approximated by a constant in the calculation of the spherical Compton signal and only causes a difference of ≤1% for clusters in the mass range [10¹⁴ M_⊙, 10¹⁵ M_⊙].

However, the hydrostatic masses used for constructing the UPP model are systematically underestimated, which means that the cluster radii are also biased. Integrating an X-ray-measured pressure profile over a biased volume leads to a biased Compton signal. We apply the rescaling methods discussed in Section 2.2 to the GNFW fits to scaled pressure profiles and correct the X-ray-measured radii of every REXCESS cluster, and correct the bias in the Compton parameter derived from X-ray measurements. We also calculate Y_sph analytically by integrating our DPP over the cluster within the radius r = xR_500c:

$$\begin{eqnarray}\begin{array}{rcl}{Y}_{\mathrm{sph}}({{xR}}_{500{\rm{c}}}) & = & \displaystyle \frac{4\pi {\sigma }_{{\rm{T}}}}{3{m}_{{\rm{e}}}{c}^{2}}{R}_{500{\rm{c}}}^{3}{P}_{500{\rm{c}}}\\ & & \times \,{\displaystyle \int }_{0}^{x}3{\left({x}^{{\prime} }\right)}^{2}{\mathbb{P}}({x}^{{\prime} }){\left[\displaystyle \frac{{M}_{500{\rm{c}}}}{3\times {10}^{14}{h}_{70}^{-1}{M}_{\odot }}\right]}^{{\alpha }_{{\rm{P}}}(x)}{{dx}}^{{\prime} },\end{array}\end{eqnarray} \tag{ 25 }$$

then we simplify the integral, getting

$$\begin{eqnarray}&&h{\left(z\right)}^{-2/3}{Y}_{\mathrm{sph}}({{xR}}_{500{\rm{c}}})=C(x){\left[\displaystyle \frac{{M}_{500{\rm{c}}}}{3\times {10}^{14}{h}_{70}^{-1}{M}_{\odot }}\right]}^{\alpha },\end{eqnarray} \tag{ 26 }$$

where α = 5/3 given by ${P}_{500{\rm{c}}}{R}_{500{\rm{c}}}^{3}$, since α_P(x) is set to be 0 in the DPP and has no contribution to α, and

$$\begin{eqnarray}\begin{array}{rcl}C(x) & = & 2.925\times {10}^{-5}I(x){h}_{70}^{-1}\,{\mathrm{Mpc}}^{2},\\ I(x) & = & {\displaystyle \int }_{0}^{x}3{\mathbb{P}}({x}^{{\prime} }){\left({x}^{{\prime} }\right)}^{2}{{dx}}^{{\prime} }.\end{array}\end{eqnarray} \tag{ 27 }$$

We use the values for the parameters P₀, c₅₀₀, α, β, γ of P_m reported in Table 2 to get I(1) = 0.554. Rewriting C(1) to the logarithmic form 10^A , we get A = −4.790, along with α = 5/3, as previously discussed for the Y_sph–M relation, which agrees well with the direct fit to the REXCESS sample after correcting for hydrostatic bias: α = 1.673 ± 0.014 and A = −4.786 ± 0.004.

In Figure 4, we plot fits for the integrated Compton signal versus cluster mass after correction for hydrostatic bias and analytical calculated Y–M relation based on DPP. For comparison, we also plot the Y–M relation reported in Arnaud et al. (2010).

Zoom In Zoom Out Reset image size

The corrected Y–M relation leads to smaller Y values at a given M compared to the UPP model, which indicates that our Y–M relation predicts a higher mass for the observed cluster given the same measured Compton signal. The new fit also shows a negligible difference from analytical results based on the DPP. The value of the best-fit slope is close to the self-similar scaling with a tiny deviation Δα = 1.673 − 5/3 = 0.006.

We find no evidence for a power-law index of the Y_sph–M_500c relation that deviates from the predictions of self-similarity, which is consistent with a small deviation from standard self-similarity of the Y_sph–mass scaling relation in Gupta et al. (2017). The disappearance of the deviation from self-similarity is mainly due to the dependence of hydrostatic bias on cluster mass. We find changes in the spherical Compton signal of clusters in the REXCESS sample after adjusting for hydrostatic bias are much less significant, <2% compared to the correction of cluster masses, which means the shift in cluster masses is the key factor for the modification on the power-law index of the Y–M relation. Notice that the relation of ${M}_{500{\rm{c}}}^{\mathrm{True}}$ versus ${M}_{500{\rm{c}}}^{{\rm{X}} \mbox{-} \mathrm{ray}}$ for the REXCESS sample we derived yields B_M ∝ M^≃0.08, equal to the shift of cluster mass after correction. To a great extent, this explains the variation—around 0.12—of the power-law index of the Y–M relation after adjusting for hydrostatic bias.

The slope of the Y–M relation being consistent with the standard self-similar model after adjusting for hydrostatic bias indicates that the studies that claim deviations from self-similarity in mass scaling relation like L_X − M and Y_X − M (e.g., Allen et al. 2003; Arnaud et al. 2007) may need to be revised as their results are also affected by similar X-ray mass biases. Additionally, the existence of self-similarity in the mass independence of our pressure model and the Y–M relation makes us more confident in extrapolating out the DPP model to higher redshifts in the calculation of the thermal SZ power spectrum by assuming a self-similarity in the redshift dependence according to the standard self-similar model. We also note that the REXCESS sample has a limited mass range. Further confirmation of self-similarity in the Y–M relation requires joint work of simulations and further observations with an extended mass range.

The tSZ power spectrum is a powerful probe of cosmology and can provide promising constraints on cosmological parameters: ${C}_{\ell }\propto {\sigma }_{8}^{7-9}$ (e.g., Komatsu & Seljak 2002; Shaw et al. 2010; Trac et al. 2011). Since clusters are the dominant source of tSZ anisotropies, due to the number density of clusters and the gas thermal pressure profile, the tSZ power spectrum can be adequately modeled by an approach referred to as the halo formalism (e.g., Cole & Kaiser 1988; Komatsu & Kitayama 1999).

The tSZ angular power spectrum at a multipole moment, $\ell$, for the one-halo term is given by

$$\begin{eqnarray}&&{C}_{\ell }^{\mathrm{tSZ}}={f}^{2}(\nu ){\int }_{z}\displaystyle \frac{{dV}}{{dz}}{\int }_{M}\displaystyle \frac{{dn}(M,z)}{{dM}}| {\tilde{y}}_{\ell }(M,z){| }^{2}{dMdz},\end{eqnarray} \tag{ 28 }$$

where $f(\nu )=x\coth (x/2)-4$ is the spectral dependence with x = h ν/(k_B T_CMB). Integrations over redshift and mass are carried out from z = 0.0 to z = 6.0 and from M = 10¹⁰ M_⊙ to M = 10¹⁶ M_⊙ respectively. For the differential halo mass function ${dn}(M,z)/{dM}$, we adopt the fitting function from Tinker et al. (2008) based on N-body simulations.

Following Komatsu & Seljak (2002), the 2D Fourier transform of the projected Compton y-parameter, ${\tilde{y}}_{l}(M,z)$ is given by

$$\begin{eqnarray}&&{\tilde{y}}_{\ell }(M,z)=\displaystyle \frac{4\pi {r}_{{\rm{s}}}}{{\ell }_{{\rm{s}}}^{2}}\displaystyle \frac{{\sigma }_{{\rm{T}}}}{{m}_{{\rm{e}}}{c}^{2}}\int {x}^{2}{P}_{{\rm{e}}}(x)\displaystyle \frac{\sin (\ell x/{\ell }_{{\rm{s}}})}{\ell x/{\ell }_{{\rm{s}}}}{dx},\end{eqnarray} \tag{ 29 }$$

with the Limber approximation (Limber 1953), which is used to relate the angular correlation function to the corresponding three-dimensional spatial clustering in an approximate way and to avoid spherical Bessel function calculations, where x = r/r_s is a scaled dimensionless radius, r_s is characteristic radius for a NFW profile defined by R_500c/c_500c, and we use average halo concentrations, c_500c, calibrated as a function of cluster mass and redshift from Diemer & Kravtsov (2015). The corresponding angular wavenumber ${\ell }_{{\rm{s}}}={d}_{{\rm{A}}}/{r}_{{\rm{s}}}$, where d_A(z) is the proper angular-diameter distance at redshift z. P_e(x) is the electron pressure we discussed in Section 2.1. The integral is carried out within a spherical region with radius R ∼ 4R_500c .

In Figure 5, we compare the measured tSZ power spectrum to the one-halo term predicted by the UPP and DPP models. Our predictions for the tSZ spectrum are made by assuming the fiducial parameters Ω_m = 0.3, Ω_Λ = 0.7, Ω_b = 0.045, h = 0.7, n_s = 0.96, σ₈ = 0.8. We use measurements of the tSZ power spectrum from the analysis of ACT (Dunkley et al. 2013), SPT (George et al. 2015), and Planck (Planck Collaboration et al. 2016b), all rescaled to 146 GHz, at which f²(ν) = 1, for direct comparison. The uncertainties in the Planck 2015 data points account for statistical and systematic errors, as well as modeling uncertainties associated with correcting for foreground contamination.

Zoom In Zoom Out Reset image size

The tSZ power spectrum derived from the original UPP model predicts much higher values than observational data while our DPP model leads to the tSZ power spectrum matching the tSZ data of Planck within 1σ uncertainty for multipoles $100\leqslant \ell \leqslant 1300$. However, the tSZ power spectrum calculated with our DPP model still shows a significant tension with ACT and SPT data at $\ell =3000$. Our work shows adjusting ICM pressure profiles for hydrostatic bias due to non-thermal pressure has a significant effect on lowering the amplitude of the power spectrum by 30%–40%. This is in agreement with other work studying the change in the shape of the tSZ power spectrum after including the effect of non-thermal pressure (e.g., Battaglia et al. 2010, 2012a; Shaw et al. 2010; Trac et al. 2011).

In the analytical calculation of the tSZ power spectrum, we extrapolate our pressure model to redshifts as high as z = 6.0, even though our pressure model is built on simulation data from a low-redshift snapshot (z = 0.1). However, we show in Figure 5 that galaxy clusters of redshift z ≥ 1.0 will not significantly affect our calculation of the tSZ power spectrum at $\ell \leqslant 1300$. The tSZ power spectrum at $\ell \geqslant 3000$ shows it is more sensitive to the high-redshift clusters, which may indicate that redshift dependence could lower the tension between our calculation of the tSZ power spectrum and high-multipole observations.