Planck/SDSS Cluster Mass and Gas Scaling Relations for a Volume-Complete redMaPPer Sample

TL;DR

This study uses Planck and SDSS data to analyze SZ gas pressure profiles and derive cluster mass and gas scaling relations for a large, volume-complete sample of optically selected clusters, revealing new insights into cluster properties.

Contribution

It provides the first detailed SZ-based pressure profiles and mass-richness relations for a large, volume-complete sample of optically selected clusters, extending previous work to lower masses.

Findings

SZ pressure profiles are well detected over a wide mass range.

SZ-based masses are ~24% lower than weak lensing masses, indicating hydrostatic bias.

Derived Y_500-M_500 relation has a slope of 1.72 +/- 0.07, consistent with previous SZ studies.

Abstract

Using Planck satellite data, we construct SZ gas pressure profiles for a large, volume-complete sample of optically selected clusters. We have defined a sample of over 8,000 redMaPPer clusters from the Sloan Digital Sky Survey (SDSS), within the volume-complete redshift region 0.100 < z < 0.325, for which we construct Sunyaev-Zel'dovich (SZ) effect maps by stacking Planck data over the full range of richness. Dividing the sample into richness bins we simultaneously solve for the mean cluster mass in each bin together with the corresponding radial pressure profile parameters, employing an MCMC analysis. These profiles are well detected over a much wider range of cluster mass and radius than previous work, showing a clear trend towards larger break radius with increasing cluster mass. Our SZ-based masses fall ~24% below the mass-richness relations from weak lensing, in a similar fashion as the "hydrostatic bias" related with X-ray derived masses. We correct for this bias to derive an optimal mass-richness relation finding a slope 1.22 +/- 0.04 and a pivot mass log(M_500/M_0)= 14.432 +/- 0.041, evaluated at a richness lambda=60. Finally, we derive a tight Y_500-M_500 relation over a wide range of cluster mass, with a power law slope equal to 1.72 +/- 0.07, that agrees well with the independent slope obtained by the Planck team with an SZ-selected cluster sample, but extends to lower masses with higher precision.

Tables

Table 1. Table 1: The recovered values of the GNFW universal pressure profile parameters, together with the best fit values.

Parameter	Mean value	Best fit
$P_0$	$6.48\pm 3.00$	$5.42$
$c_{500}$	$1.65\pm 0.95$	$1.12$
$\alpha$	$0.87\pm 0.25$	$0.81$
$\beta$	$3.19\pm 0.66$	$3.47$

Table 2. Table 2: Richness range, number of clusters, mean richness and mean masses recovered from the SZ signal of the cluster subsamples studied.

Richness range	$N$		$\left[{10}^{14}{M}_{\odot }\right]$
$[136,\mathrm{\hspace{0.17em}200})$	19	157.6	$6.48\pm 0.66$
$[93,\mathrm{\hspace{0.17em}136})$	68	109.6	$4.52\pm 0.42$
$[63,\mathrm{\hspace{0.17em}93})$	293	74.3	$2.77\pm 0.25$
$[43,\mathrm{\hspace{0.17em}63})$	902	50.8	$1.84\pm 0.18$
$[29,\mathrm{\hspace{0.17em}43})$	2308	34.6	$1.04\pm 0.13$
$[20,\mathrm{\hspace{0.17em}29})$	4440	23.8	$0.54\pm 0.12$

Table 3. Table 3: The confidence values of the mass–richness relation parameters and of the bias factor ( 1 − b ) 1 𝑏 (1-b) , together with the best fit. α M | λ subscript 𝛼 conditional 𝑀 𝜆 \alpha_{M|\lambda} is the slope of the mass–richness relation, and log 10 ⁡ M 0 subscript 10 subscript 𝑀 0 \log_{10}M_{0} is the pivot mass, evaluated at λ 0 = 60 subscript 𝜆 0 60 \lambda_{0}=60 and given in terms of M 500 ​ c subscript 𝑀 500 𝑐 M_{500c} . The stacking analysis performed is not able to constrain the value of σ M | λ subscript 𝜎 conditional 𝑀 𝜆 \sigma_{M|\lambda} , which is largely determined by the prior, and thus not shown.

Parameter	Mean value	Best fit
${\log }_{10}\left(M_0/M_{\odot }\right)$	$14.435\pm 0.040$	$14.446$
${\alpha }_{M|\lambda }$	$1.22\pm 0.04$	$1.21$
$\left(1-b\right)$	$0.76\pm 0.05$	$0.75$
${\sigma }_{M|\lambda }$	-	-

Equations

$$\lambda =\sum p_{\mathrm{mem}}{\theta }_L{\theta }_R,$$

$$\frac{\Delta T_{\mathrm{SZ}}}{T_{\mathrm{CMB}}}=g(x)y,$$

$$g(x)=\left(x\frac{e^x+1}{e^x-1}-4\right).$$

$$y=\frac{{\sigma }_T}{m_e{c}^2}{\int }_0^{\infty }P(l)\mathrm{d}l,$$

$$Y={\int }_{\Omega }y\mathrm{d}\Omega =D_A^{-2}\frac{{\sigma }_T}{m_e{c}^2}{\int }_0^{\infty }\mathrm{d}l{\int }_{A_{\mathrm{clu}}}P(l)\mathrm{d}A,$$

$$y(r)=\frac{{\sigma }_T}{m_e{c}^2}{\int }_{-\infty }^{\infty }P\left(\sqrt{{r^{\prime }}^2+r^2}\right)\mathrm{d}{r}^{\prime },$$

$$Y(R)={\int }_0^{R}2\pi y(r)r\mathrm{d}r,$$

$$Y^{\mathrm{sph}}(R)=\frac{{\sigma }_T}{m_e{c}^2}{\int }_0^{R}4\pi P(r)r^2\mathrm{d}r.$$

$$Y[{\mathrm{arcmin}}^2]=D_A(z{)}^{-2}{\left(\frac{60\times 180}{\pi }\right)}^{2}Y[{\mathrm{Mpc}}^2]$$

$$\mathbb{P}(x)=\frac{P_0}{(c_{500}x{)}^{\gamma }{\left[1+(c_{500}x{)}^{\alpha }\right]}^{(\beta -\gamma )/\alpha }},$$

$$P(x)=P_{500}{\left(\frac{M_{500}}{3\times 10^{14}{h}_{70}^{-1}{M}_{\odot }}\right)}^{{\alpha }_p}\mathbb{P}(x),$$

$$P_{500}=$$

$$\times {\left(\frac{M_{500}}{3\times 10^{14}{h}_{70}^{-1}{M}_{\odot }}\right)}^{2/3}{\mathrm{h}}_{70}^2\mathrm{keV}{\mathrm{cm}}^{-3},$$

$$\frac{4}{3}\pi 500{\rho }_c(z)r_{500}^3=M_{500},$$

$${\rho }_c(z)=\frac{3H_0^{2}E(z{)}^2}{8\pi G},$$

$$\ln \left(f_{\mathrm{gas},500}E(z{)}^{-3/2}\right)=$$

$$+\left(0.21\pm 0.03\right)\ln \left(\frac{M_{500}}{2\times 10^{14}{M}_{\odot }}\right).$$

$$M_{\mathrm{gas},500}={\int }_0^{r_{500}}{\mu }_e{m}_u{n}_e(r)4\pi r^2\mathrm{d}r,$$

$${\log }_{10}\left(\frac{M_{500}E(z)}{h_{70}^{-1}{M}_{\odot }}\right)=$$

$$+\left(1.67_{-0.10}^{+0.14}\right){\log }_{10}\left(\frac{k_{B}T}{\mathrm{keV}}\right),$$

$$M_{\mathrm{HE},500}=(1-b)M_{500},$$

$$⟨M|\lambda ⟩=M_0{\left(\frac{\lambda }{{\lambda }_0}\right)}^{{\alpha }_{M|\lambda }},$$

$$P(M|{\lambda }^{\mathrm{obs}})=\int P(M|\lambda )P(\lambda |{\lambda }^{\mathrm{obs}})\mathrm{d}\lambda ,$$

$$P(\lambda |{\lambda }^{\mathrm{obs}})=\frac{1}{\sqrt{2\pi {\sigma }_{\ln \lambda |{\lambda }^{\mathrm{obs}}}^2}}\exp [-x^2(\lambda )],$$

$$x(\lambda )=\frac{\ln \lambda -\ln {\lambda }^{\mathrm{obs}}}{\sqrt{2{\sigma }_{\ln \lambda |{\lambda }^{\mathrm{obs}}}^2}},$$

$$⟨M⟩=\frac{1}{N_i}\sum _{j=1}^{N_i}\int MP(M|{\lambda }_j^{\mathrm{obs}})\mathrm{d}M.$$

$$\ln \mathcal{L}(\mathbf{\varphi })=\sum _{k=1}^6\ln {\mathcal{L}}^k({\mathbf{y}}^k|{\mathbf{\varphi }}^k),$$

$$\ln {\mathcal{L}}^k({\mathbf{y}}^k|{\mathbf{\varphi }}^k)\propto -\frac{1}{2}{\chi }^2({\mathbf{y}}^k,{\mathbf{\varphi }}^k,\mathcal{C}),$$

$${\chi }^2({\mathbf{y}}^k,{\mathbf{\varphi }}^k,\mathcal{C})=\left({\mathbf{y}}^k-\mathbf{\mu }({\mathbf{\varphi }}^k)\right){\mathcal{C}}^{-1}{\left({\mathbf{y}}^k-\mathbf{\mu }({\mathbf{\varphi }}^k)\right)}^T,$$

$${\mathcal{C}}_{ij}(M_{500})=\frac{1}{N_R-1}\sum _{n=1}^{N_R}(y_i^n-⟨y_i⟩)(y_j^n-⟨y_j⟩).$$

$$\ln {\mathcal{L}}_{\mathrm{bias}}=\ln {\mathcal{L}}_{S16}+\ln {\mathcal{L}}_{M16},$$

$${\chi }^2\left[(1-b)\right]=\left({\mathbf{M}}_{500}-{\mathbf{M}}_{500}^{\text{model}}\right){\mathcal{C}}^{-1}{\left({\mathbf{M}}_{500}-{\mathbf{M}}_{500}^{\text{model}}\right)}^T,$$

Full text

Our SZ-based masses fall $\sim$24% below the mass–richness relations from weak lensing, in a similar fashion as the “hydrostatic bias” related with X-ray derived masses. We correct for this bias to derive an optimal mass–richness relation finding a slope $1.22\pm 0.04$ and a pivot mass $\log_{10}(M_{500}/M_{\odot})=14.432\pm 0.041$, evaluated at a richness $\lambda=60$. Finally, we derive a tight $Y_{500}$–$M_{500}$ relation over a wide range of cluster mass, with a power law slope equal to $1.72\pm 0.07$, that agrees well with the independent slope obtained by the Planck team with an SZ-selected cluster sample, but extends to lower masses with higher precision.

††pagerange: Planck/SDSS Cluster Mass and Gas Scaling Relations for a Volume-Complete redMaPPer Sample–7

1 Introduction

Galaxy clusters are powerful cosmological probes that provide complementary constraints in the era of “Precision Cosmology”. They contribute accurate consistency checks and unique new competitive constraints because of the well understood cosmological sensitivity of their numbers and clustering (Jain et al., 2013; Huterer et al., 2015; Dodelson & Park, 2014; Pouri et al., 2014; Pan & Knox, 2015). The growth of structures has led the observational evidence to support dark energy dominance today, in combination with complementary constraints (Efstathiou et al., 1990; Lahav et al., 1991; Bahcall, 2000; Allen et al., 2011; Carroll et al., 1992; Ostriker & Steinhardt, 1995; Bahcall & Fan, 1998). To realize their full cosmological potential, large, homogeneous samples of clusters are now being constructed out to $z\simeq 1$ with weak lensing based masses, in particular the Subaru/HSC and JPAS surveys (Oguri et al., 2017; Benitez et al., 2014; Jimeno et al., 2015, 2017).

Currently the best direct lensing masses are limited to relatively small subsamples of X-ray and Sunyaev-Zel’dovich (SZ) selected clusters, totalling $\sim$100 clusters (Umetsu et al., 2014; Zitrin et al., 2015; Okabe & Smith, 2016). One of the main efforts is focused on defining scaling relations between clusters with such weak lensing masses and the more widely available X-ray, SZ and/or optical richnesses with the reasonable expectation that these relations may provide mass proxies in the absence of direct lensing masses. Such proxies have a physical basis for clusters that appear to be virialised, so that X-ray temperature and emissivity profiles can provide virial masses under hydrostatic equilibrium. Independently, the SZ distortion of the cosmic microwave background (CMB) spectrum relates the density and temperature of cluster gas through inverse Compton scattering, and hence naturally anticipated to scale approximately with cluster mass. The cleanest mass proxy is arguably provided by the number of member galaxies, the so called richness, implicit in the assumption that the dominant cluster dark matter is collisionless like galaxies, and indeed the mass and the richness are found to be closely proportional (Bahcall, 1977; Girardi et al., 2000).

The cluster mass–richness relation, crucial in any attempt to use large number of clusters detected in the optical to constrain cosmological parameters, has been estimated in the past decade using cluster catalogues derived from the Sloan Digital Sky Survey (SDSS, Gunn et al. 2006) data, like the MaxBCG (Koester et al., 2007), the GMBCG (Hao et al., 2010), the WHL12 (Wen et al., 2012), or the more sophisticated redMaPPer cluster catalogue, both in its SDSS (Rykoff et al., 2014) and DES (Rykoff et al., 2016) versions. This relation can be estimated directly obtaining cluster masses from X-rays, weak lensing, SZ effect or velocity dispersion measurements in clusters (Johnston et al., 2007; Andreon & Hurn, 2010; Saro et al., 2015; Sereno et al., 2015; Saro et al., 2016; Simet et al., 2016; Mantz et al., 2016; Melchior et al., 2016), or indirectly, using numerical simulations (Angulo et al., 2012; Campa et al., 2017) or comparing the observed abundances or clustering amplitudes with model predictions (Rykoff et al., 2012; Baxter et al., 2016; Jimeno et al., 2017).

The dynamical evolution and growth of galaxy clusters are driven by the dominant dark matter, but the relevant observables depend on the physical state of the baryons. Hence scaling relations between clusters observables and mass are not direct, but have been predicted to follow physically self-similar relations (Kaiser, 1986; Kravtsov & Borgani, 2012) that have been tested observationally and with hydrodynamical N-body simulations. Specifically, the integrated thermal SZ effect (Sunyaev & Zeldovich, 1972), the X-ray luminosity, and the temperature are predicted to scale with the mass of the galaxy clusters as $Y\propto M^{5/3}$, $L_{X}\propto M^{4/3}$ and $T\propto M^{2/3}$, respectively. While simulations agree with the self-similar model (White et al., 2002; da Silva et al., 2004; Motl et al., 2005; Nagai, 2006; Wik et al., 2008; Aghanim et al., 2009), X-ray and SZ observations have uncovered departures from self-similarity that may be explained by complications due to cluster mergers, including shocked gas, cool gas cores, and energy injection from active galactic nuclei (AGN)(Voit, 2005; Arnaud et al., 2005, 2007; Pratt et al., 2009; Vikhlinin et al., 2009; De Martino & Atrio-Barandela, 2016). Differences between the observations and purely gravitationally predicted scaling relations then provide insights into the interesting physics of the intracluster medium (Bonamente et al., 2008; Marrone et al., 2009; Arnaud et al., 2010; Melin et al., 2011; Andersson et al., 2011; Comis et al., 2011; Czakon et al., 2015).

In practice, samples of strong SZ-selected clusters that are also bright X-ray sources are currently being used to calibrate the SZ–mass relation (Arnaud et al., 2010; Planck Collaboration et al., 2013a, 2014d; Saro et al., 2015), but, since such clusters are often out of hydrostatic equilibrium for the reasons mentioned above, an SZ–mass scaling relation requires a correction for “hydrostatic mass bias” (Nagai et al., 2007b; Zhang et al., 2010; Shi & Komatsu, 2014; Sayers et al., 2016). This hot gas related bias can be broadened by other systematics like object selection process or by temperature inhomogeneities in X-ray measurements. Another approach to calibrate the masses of the cluster sample is to stack clusters in terms of richness and measure the SZ signal as a function of richness. This was done first by Planck Collaboration et al. (2011) using the MaxBCG catalogue and stacking the Planck data, and more recently by Saro et al. (2016) using an initial sample of 719 DES clusters with South Pole Telescope (SPT) SZ data and assuming various priors to extract the SZ signal.

In this work, we extract the Planck SZ signal from $\sim$ 8,000 redMaPPer clusters identified in the SDSS, that have allowed us to previously define accurate clustering and density evolution measurements in the redshift range $0.100<z<0.325$, as described in detail in Jimeno et al. (2017). Here we take this well defined cluster sample and stack the Planck multi-frequency data over a wide range of cluster richness. We only need to assume a weak prior for the global gas fraction using X-ray measurements to simultaneously derive a more “self-sufficient” method to derive SZ pressure profiles and the corresponding mean cluster masses binned by richness. Comparing masses derived this way with those expected from weak lensing mass–richness relations found in the literature, we derive the level of intrinsic bias for our sample and we then derive both a debiased mass–richness and a $Y_{500}$–$M_{500}$ relation describing our observational results.

This paper is organised as follows. In Sec. 2 we describe the data that we use in our analysis, namely the redMaPPer cluster catalogue and Planck HFI maps. We present the basic theory associated to the SZ effect in Sec. 3, as well as the mass–richness relation and the different models that we consider. In Sec. 4 we process the Planck data to obtain SZ maps given in terms of the Compton parameter $y$ and use them to constrain, through a joint likelihood analysis, the universal pressure profile parameters and the mean masses of the cluster subsamples considered. In Sec. 5 we make an estimation of the value of the mass bias and obtain the optimal mass–richness relation able to describe our bias-corrected masses. Finally, we use all the results obtained in the previous sections to derive a $Y_{500}$–$M_{500}$ relation in Sec. 6, and present our conclusions in Sec. 7.

Throughout this paper we adopt a fiducial flat $\Lambda{\rm CDM}$ cosmology with a matter density $\Omega_{m}=0.3$ and a Hubble parameter with a value today of $H_{0}=70\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$. We also consider ${h\equiv H_{0}/(100\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1})}$ and ${h_{70}\equiv H_{0}/(70\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1})}$. Cluster masses are given in terms of an spherical overdensity $\Delta_{c(m)}$ with respect to the critical (mean) density of the Universe, $M_{\Delta{c(m)}}=(4/3)\pi\Delta\rho_{c(m)}(z)r_{\Delta}^{3}$. Unless stated otherwise, we refer to $M_{500c}$ as $M_{500}$.

2 Data

2.1 redMaPPer cluster catalogue

The pioneering Sloan Digital Sky Survey (SDSS, Gunn et al. 2006) combines photometric and spectroscopic observations and has mapped the largest volume of the Universe in the optical to date, covering around 14,000 deg2 of the sky. The information obtained has been made public periodically via Data Releases (DR), and many different cluster catalogues have been carefully constructed using very different cluster-finder algorithms (Koester et al., 2007; Hao et al., 2010; Wen et al., 2012). We focus our analysis on one of these SDSS-based catalogues, the red-sequence Matched-filter Probabilistic Percolation cluster catalogue (redMaPPer, Rykoff et al., 2014), based on SDSS DR8 photometric data, as given in the public 6.3 version (Rykoff et al., 2016). We use this catalogue because it offers a very low rate of projection effects ($<5\%$, according to the thorough analysis of Rykoff et al. (2014)). We have previously compared this catalogue with other cluster catalogues and found the redMaPPer catalogue to provide highly consistent results when comparing gravitational redshift and the level of magnification bias (Jimeno et al., 2015) and more recently we have also constructed a spectroscopically sample of redMaPPer clusters to construct the most precise cluster correlation functions and mass distributions to date (Jimeno et al., 2017), finding a high degree of consistency in relation to $\Lambda$CDM based predictions from the forefront MXXL simulations (Angulo et al., 2012).

The redMaPPer cluster finder algorithm, based on Rozo et al. (2009b) and Rykoff et al. (2012), is able to find potential clusters within the raw photometric data, and provides an estimation of the sky position of the five most probable central galaxies (CGs), the redshift $z_{\lambda}$, and the richness $\lambda$ of each cluster. In this catalogue, the richness estimations relies on a multi-colour self-training procedure that calibrates the red-sequence as a function of redshift. The richness of a cluster is defined as:

[TABLE]

where $p_{\rm mem}$ is the membership probability of each galaxy found near the cluster, and $\theta_{L}$ and $\theta_{R}$ are luminosity and radius-dependent weights. A more in-depth explanation of the algorithm features can be found in Rykoff et al. (2014) and Rozo et al. (2015b).

The resulting redMaPPer catalogue covers an effective area of 10,401 deg2, and contains 26,111 clusters in the $0.08\leqslant z_{\rm photo}\leqslant 0.55$ redshift range. Finally, it should be mentioned that this catalogue is volume-complete up to $z\lesssim 0.33$, and has a richness cutoff of $\lambda/S(z)>20$, where $S(z)$ is the “scale factor” that relates the richness with the number of observed galaxies above the magnitude limit of the survey, $\lambda/S(z)$. This detection threshold corresponds to a mass limit of approximately $M_{500}\geqslant 0.7\times 10^{14}\,M_{\odot}$.

2.2 Planck SZ data

We combine the above optically selected cluster sample from SDSS with the all-sky temperature maps derived by the Planck space mission (Planck Collaboration et al., 2014a). Although these maps have already been used to construct catalogues of SZ sources (Planck Collaboration et al., 2014b) and an all-sky Compton $y$ parameter map (Planck Collaboration et al., 2014e), we reprocess them for our own purposes.

To obtain the Compton parameter maps required in our analysis, we use the Planck full mission High-Frequency Instrument maps (HFI, Planck Collaboration et al. 2014c) at 100, 143, 217, and 353 GHz. These maps are provided in HEALPix format (Górski et al., 2005), with a pixelisation of $N_{\rm side}=$ 2,048, which correspond to a pixel resolution of $\sim 1.7$ $\mathrm{arcmin}$. The Planck effective beams for each of the 100, 143, 217, and 353 GHz channels can be approximated by circular Gaussians with FWHM values of 9.66, 7.27, 5.01 and 4.86 $\mathrm{arcmin}$, respectively. To compute the contribution of the SZ signal in the Planck temperature maps, we also make use of the spectral transmission information of each of these frequency channels, as given in Planck Collaboration et al. (2014c).

3 Model

3.1 Thermal Sunyaev-Zel’dovich effect

Here we briefly introduce the equations that describe the thermal Sunyaev-Zel’dovich (SZ) effect. For a derivation, we refer the reader to the papers of Sunyaev & Zel’dovich (1980); Rephaeli (1995), or the more recent work by Birkinshaw (1999), Carlstrom et al. (2002) and Diego et al. (2002). Ignoring relativistic corrections, the SZ spectral distortion of the CMB, expressed as temperature change, is:

[TABLE]

where $x=(h\,\nu)/(k_{B}\,T_{\rm CMB})$ is the dimensionless frequency, $y$ is the Compton parameter, and

[TABLE]

The Compton parameter $y$ is equal to the optical depth, $\tau_{e}$, times the fractional energy gain per scattering, and is given by:

[TABLE]

where $\sigma_{T}$ is the Thomson cross section and $P$ is the intracluster pressure produced by free electrons. Integrating over the solid angle of the cluster one obtains the integrated Compton parameter:

[TABLE]

where $A_{\rm clu}$ is the area of the cluster in the plane of the sky, and $D_{A}(z)$ is the angular diameter distance at redshift $z$.

Assuming an spherical model for the cluster, we have that the Compton parameter $y$ at a distance $r$ from the center of the cluster is equal to:

[TABLE]

and thus the integrated Compton parameter $Y$, obtained integrating to a distance $R$ from the center of the cluster, is given by:

[TABLE]

which has units of $\mathrm{Mpc}^{2}$. It should be noted that, as $y$ is a projected along the line of sight quantity, $Y$ is the so called “cylindrical” integrated Compton parameter $Y^{\rm cyl}$, and not the “spherical” integrated Compton parameter, which would be obtained directly from the pressure profile doing:

[TABLE]

In practice, we work with $Y^{\rm cyl}$ when dealing with observations, as this is the quantity that can be measured from the data, and we use $Y^{\rm sph}$ when dealing with models. Once a pressure profile has been adopted, any measurement of $Y^{\rm cyl}(n\,r_{500})$ can be straightforwardly converted in terms of $Y^{\rm cyl}(r_{500})$, and the latter to $Y^{\rm sph}(r_{500})$. We refer the reader to appendix A of Melin et al. (2011) for a more detailed explanation of how to convert between definitions.

Finally, as $y$ is dimensionless, $Y$ can also be expressed in units of $\mathrm{arcmin}^{2}$:

[TABLE]

From now on, we refer to $Y^{\rm sph}(r_{500})$ as $Y_{500}$, given in $\mathrm{Mpc}^{2}$ units.

3.2 Pressure profile

In this work we adopt the generalised NFW (GNFW) “universal pressure profile” proposed by Nagai et al. (2007b), that has a flexible double power-law form:

[TABLE]

where $x=r/r_{500}$ is the scaled dimensionless physical radius. The physical pressure is given by:

[TABLE]

where:

[TABLE]

and ${\alpha_{p}=0.12}$ accounts for the deviation from the self-similar scaling model (Arnaud et al., 2010). The value of $r_{500}$ is given by:

[TABLE]

where $\rho_{c}(z)$ is the critical density of the Universe at redshift $z$, defined as:

[TABLE]

and $E(z)^{2}=\Omega_{m}(1+z)^{3}+(1-\Omega_{m})$.

From Eq. 10, it is clear that the slopes of the pressure profile are given, at different $r_{500}$-scaled distances, by $\gamma$ for $x\ll 1/c_{500}$, $\alpha$ for $x\sim 1/c_{500}$, and $\beta$ for $x\gg 1/c_{500}$. In our analysis and following the approach by Planck Collaboration et al. (2013a), we leave $P_{0}$, $c_{500}$, $\alpha$, and $\beta$ as free parameters. The low resolution of the Planck data does not have the power to constrain $\gamma$, so we fix it to $\gamma=0.31$, value obtained by Arnaud et al. (2010) from a sample of 33 XMM-Newton X-ray local clusters in the $r<r_{500}$ range.

3.3 Gas fraction

To improve our analysis, we use established results regarding the global gas fraction $f_{\mathrm{gas}}$ in clusters, particularly, those by Pratt et al. (2009), who derived a mass–gas fraction relation using precise hydrostatic mass measurements of 41 Chandra and XMM-Newton clusters (Vikhlinin et al., 2006; Arnaud et al., 2007; Sun et al., 2009), which is also in good agreement with the results obtained from the REXCESS sample (Böhringer et al., 2007). According to their analysis, these clusters, whose masses range from $10^{13}\,M_{\odot}$ to $10^{15}\,M_{\odot}$, follow the mean mass–gas fraction relation:

[TABLE]

To compute the gas fraction we first need to compute the gas mass:

[TABLE]

where $\mu_{e}=1.15$ is the mean molecular weight per free electron, $m_{u}$ is the atomic mass unit, and $n_{e}(r)$ is the electron density. Because the intra-cluster pressure is given by $P(r)=n_{e}(r)\,k_{B}\,T$, assuming an isothermal model for the cluster one can directly derive $M_{{\rm gas},500}$ from the adopted pressure profile (Eq. 11).

For the temperature, we use the mean mass–temperature relation given by Lieu et al. (2016):

[TABLE]

which was obtained combining weak lensing mass estimates with Chandra and XMM-Newton temperature data of 38 clusters from the XXL survey (Pacaud et al., 2016), 10 clusters from the COSMOS survey (Kettula et al., 2013), and 48 from the Canadian Cluster Comparison Project (CCCP, Mahdavi et al. 2013; Hoekstra et al. 2015), spanning a temperature range $T\simeq 1-10$ keV.

It is worth mentioning that if an isothermal model is assumed and we consider that $M_{\rm gas}\propto f_{b}\,M$, where $M$ is the total cluster mass and $f_{b}$ is the baryon gas fraction, from Eq. 5 we have that the integrated Compton parameter scales as $Y\propto f_{b}\,M\,T\,D_{A}^{-2}$. However, even clusters in hydrostatic equilibrium are not strictly isothermal, and temperatures are commonly observed to drop by a factor of $\sim$2 below a radius of $r\lesssim 100-200$ kpc because of strong radiative cooling, described best by a broken power law with a transition region (Vikhlinin et al., 2006). In any case, these scales are not resolved by Planck and in our analysis the assumption that the temperature is constant is a good approximation for the radial scales considered in this work.

3.4 Miscentering

In the redMaPPer catalogue, for each cluster the 5 most probable central galaxies (CGs) are provided with their corresponding centering probabilities. Usually, there is one CG with a much higher probability of being the real CG than the other 4, so we consider the most probable CG to be the center of the cluster. In any case, it is now known that, because clusters are still evolving systems, CGs do not always reside at the deepest part of the DM halo potential well (von der Linden et al., 2007), but sometimes have high peculiar velocities, are displaced with respect to the peak of the X-ray emission (Rozo & Rykoff, 2014), or are wrongly identified satellite galaxies (Skibba et al., 2011).

In stacked measurements on clusters, miscentering is one of the main sources of noise, and should be taken into account. When modelling the SZ signal coming from stacked samples of clusters, we introduce this effect considering the results obtained by Johnston et al. (2007), who found a CG-center offset distribution that could be fitted by a 2D Gaussian with a standard deviation of $\sigma=0.42\,h^{-1}\,\mathrm{Mpc}$ for the CGs that were not accurately centered, that occurs between 20 and 40 per cent of the time as a function of cluster richness, with a probability $p_{mc}(\lambda)=(2.13+0.046\,\lambda)^{-1}$.

However, it should be noted that this value of $0.42\,h^{-1}\,\mathrm{Mpc}$ is about 2 $\mathrm{arcmin}$ at $z\sim 0.2$, scale well below the resolution of the Planck data we work with, so we do not expect this miscentering to introduce a high level of noise in our stacked measurements of the SZ effect.

3.5 Mass bias

Usually referred to as hydrostatic equilibrium (HE) masses, in their derivation there is an implicit assumption that the pressure is purely thermal. However, we may expect a non-negligible contribution to the total pressure from bulk and turbulent gas motions related to structure formation history, magnetic fields, and AGN feedback (Shi & Komatsu, 2014; Planck Collaboration et al., 2014d). Such non-thermal contributions to the total pressure would therefore cause masses estimated using X-ray or SZ observations to be biased low with uncertain estimates ranging between 5% to 20% (Nelson et al., 2014; Nagai et al., 2007a; Rasia et al., 2006; Sembolini et al., 2013).

We simply relate the HE mass estimates $M_{{\rm HE,}500}$ obtained from our SZ observations to true masses $M_{500}$ through a simple mass independent bias:

[TABLE]

where $(1-b)$ is the bias factor. This term can include not only the bias coming from departures from HE, but from observational systematics or sample selection effects.

3.6 Mass–richness relation

In order to explore the connection between the mass and the optical richness $\lambda$ in clusters, i.e., the number of galaxies contained within them, one needs to assume a form to describe the relation between cluster richness and mass. We consider the standard power law cluster mass–richness mean relation:

[TABLE]

where $M_{0}$ is a reference mass at a given value of $\lambda=\lambda_{0}$, and $\alpha_{M|\lambda}$ is the slope of the mass–richness relation. In our case, we consider $\lambda_{0}=60$.

To compute the mean masses of our cluster subsamples, we first consider the probability $P(M|\lambda^{\rm obs})$ of having a given value of the mass $M$ for a cluster with $\lambda^{\rm obs}$:

[TABLE]

with $P(M|\lambda)$ a delta function, as the relation between mass and richness is given by Eq. 19. Following the usual approach (Lima & Hu, 2005), we consider that $P(\lambda|\lambda^{\rm obs})$ follows a log-normal distribution:

[TABLE]

where:

[TABLE]

and $\sigma_{\ln\lambda|\lambda^{\rm obs}}$ is the fractional scatter on the halo richness at fixed observed richness, which is assumed to be constant with cluster redshift and richness. Because $P(M|\lambda)$ is a delta function, and considering Eq. 19, we also have that $\sigma^{2}_{\ln\lambda|\lambda^{\rm obs}}=\sigma^{2}_{\ln M|\lambda^{\rm obs}}$. So, the mean mass $\left<M\right>$ of one of the richness bins considered, with $\lambda\in[\lambda^{\rm obs}_{i},\,\lambda^{\rm obs}_{i+1}]$ and containing $N_{i}$ clusters, is given by:

[TABLE]

Finally, when dealing with mean values we can work in terms of $\left<\ln M|\lambda\right>$ instead of $\left<M|\lambda\right>$, with $\left<\ln M|\lambda\right>=\ln\left<M|\lambda\right>-0.5\,\sigma^{2}_{\ln M|\lambda^{\rm obs}}$. We refer the reader to Rozo et al. (2009a) and Simet et al. (2016) for a discussion of this transformation. We choose our parametrisation (Eq. 19) because the resulting mean relation is less affected by the uncertainty in $\sigma_{\ln M|\lambda^{\rm obs}}$.

Since it was made public, there have been multiple attempts to constrain in different ways the parameters of this relation using the redMaPPer cluster catalogue (Rykoff et al., 2012; Baxter et al., 2016; Farahi et al., 2016; Li et al., 2016; Melchior et al., 2016; Miyatake et al., 2016; Saro et al., 2015, 2016; Simet et al., 2016; Jimeno et al., 2017). Although some of these works introduced a redshift dependence in the mass–richness relation, it was weakly constrained in all cases, and compatible with no redshift evolution at all. Given the small redshift range in which we work, redshift evolution is not important for our analysis and we refer our result to the mean redshift of our sample, $z=0.245$.

4 Pressure profiles and mass estimation

4.1 Planck data processing

We first divide the redMaPPer cluster catalogue in 6 independent log-spaced richness bins, and take all those clusters that reside within the $0.100<z<0.325$ volume-complete redshift region. This leaves a total of 8,030 clusters, distributed in number and mean richness as shown in Table 4.3.

Then, for each cluster subsample, we produce and stack the $\nu=$ 100, 143, 217 and 353 GHz 2.5 deg $\times$ 2.5 deg Planck maps associated to the clusters in each subsample and produce the corresponding SZ maps following a technique similar to the one used in Planck Collaboration et al. (2016), based on internal linear combinations (ILC) of the four different HFI maps. In our case we do not use the 70 GHz SZ map, as we prefer to smooth all the maps to a common higher 10 $\mathrm{arcmin}$ resolution instead. We also use the $M_{353}-M_{143}$ combination to clean the $M_{217}$ map, where $M_{\nu}$ is the Planck map at frequency $\nu$. To convert from $\Delta T_{\rm SZ}$ to $y$ units (Eq. 2), we compute the different effective spectral responses integrating the expected SZ spectrum (Eq. 3) over each Planck bandpass.

Our final SZ maps of the cluster subsamples considered, given in terms of the Compton parameter $y$ and shown in Fig. 1, are produced as a combination of the SZ100 and SZ143 maps, weighting them by the inverse of the variance of each map. This particular combination has been proposed by the Planck team to maximise the signal-to-noise of the SZ component whilst minimising the contamination from Galactic emission and extragalactic infrared emission within clusters (Planck Collaboration et al., 2016).

4.2 Likelihood analysis

We now combine the data from Planck and the constraints imposed by the gas fractions to perform a likelihood analysis that enables us to constrain the pressure profile parameters and the mean masses of the 6 cluster subsamples. We explore the values of our 4+6 dimensional model $\bm{\phi}=(P_{0},\,c_{500},\,\alpha\,,\beta,\,M_{500}^{1},\,...\,,\,M_{500}^{6})$ through a Monte Carlo Markov Chain (MCMC) analysis.

For a given value of $r_{500}(M_{500},\,z)$ we measure the Compton parameter $y$ within a disk of radius $x=0.35$ and in 6 annulus given by radii $x_{i}$ and $x_{i+1}$, where the $x_{i}=r_{i}/r_{500}$ values are log-spaced between 0.35 and 3.5. This results in a $y$ vector of 7 values. Then, to account for the background we subtract the mean value of the signal obtained from an annulus of radii $x_{\rm out}$ and $x_{\rm out}+x_{\rm FWHM}$, where $x_{\rm out}=3.5$, and $x_{\rm FWHM}=\theta_{\rm FWHM}\,D_{A}(z)/r_{500}$ corresponds in $x$-space to the $\theta_{\rm FWHM}=10$ $\mathrm{arcmin}$ FWHM effective resolution of the SZ maps. The values of $r_{500}$ used to both model the signal and measure it from the data are obtained from $M_{500}$ through Eq. 13.

The log-likelihood employed has the form:

[TABLE]

where $\bm{\phi}^{k}=(P_{0},\,c_{500},\,\alpha\,,\beta,\,M_{500}^{k})$, $\bm{y}^{k}$ is the data vector obtained from the cluster subsample $k$, and

[TABLE]

where $\mathcal{C}$ is the covariance matrix, and

[TABLE]

with $\bm{\mu}$ the model values drawn from $\bm{\phi}^{k}$. To model the signal, for each $\bm{\phi}^{k}$ configuration we produce mock maps of the Compton parameter $y$ as a function of redshift. Then, we mimic miscentering effects adding to the mock maps the same maps smoothed with a 2D Gaussian of width $\sigma_{mc}=0.42\,h^{-1}\,\mathrm{Mpc}$, and weighted by $p_{mc}(\lambda)=(2.13+0.046\,\lambda)^{-1}$. Finally, we produce a weighted map integrating over the redshift distribution of the subsample considered, convolve it with a 10 $\mathrm{arcmin}$ FWHM Gaussian, and perform the same measurements made in the Planck data maps.

It should be noted that the mock $y$ maps that we create to fit the observed signal are generated from the pressure profile as given by Eq. 11 and for a given total model mass $M^{k}_{500}$. Hence, to compute $\ln\mathcal{L}^{k}(\bm{y}^{k}|\bm{\phi}^{k})$ we do not rely on the gas mass or $f_{gas}$.

The covariance matrixes $\mathcal{C}$ are estimated from $N_{R}=$ 1,000 patches randomly chosen within the redMaPPer footprint, where the same measurement described above is done. As this measurement depends on the $M_{500}$ value proposed, the $7\times 7$ covariance matrix is recomputed each time as:

[TABLE]

As an example, the covariance matrix obtained considering a mass of $M_{500}=5\times 10^{14}M_{\odot}$ is shown in Fig. 2.

Finally, we consider the gas fraction constraints introducing a Gaussian prior $\sim\mathcal{N}(f_{\rm gas},\,\sigma_{f_{\rm gas}})$, where $f_{\rm gas}$ is estimated from the results of Pratt et al. (2009), as explained above. In each MCMC step we compute, following the procedure described in Sec. 3.3, the 6 gas fractions associated to a given set $\bm{\phi}$ of pressure profile parameters and masses, and then use Eq. 15 to model the expected gas fraction $f_{\rm gas}$ for each value of $M_{500}^{k}$, which we use for the prior. To estimate $\sigma_{f_{\rm gas}}$, we add in quadrature the errors derived from the uncertainties on both the mass–gas fraction (Eq. 15) and the mass–temperature (Eq. 17) relations, which are obtained propagating through a Monte Carlo (MC) method. Because we have decided to be as conservative as possible on the relations employed and the resulting uncertainties in this analysis are large, we notice that the contribution that this prior has in the final estimated values of the pressure profile is small, only limiting those models where the gas fraction $f_{{\rm gas},500}$ takes values below 0.05 or above 0.20 for masses in the $\sim 10^{14}-10^{15}\,M_{\odot}$ range.

For all $P_{0}$, $c_{500}$, $\alpha$, $\beta$, $M_{500}^{1},\,...\,,\,M_{500}^{6}$, we consider flat uninformative priors $\sim\mathcal{U}(\infty,\,-\infty)$, allowing for a wide range of different model-masses configurations.

4.3 Results

The derived posterior probabilities of the GNFW universal pressure profile parameters $P_{0}$, $c_{500}$, $\alpha$ and $\beta$ are displayed in Fig. 3. To compute the center (mean) and the scale (dispersion) of the marginalised posterior distributions, we use the robust estimators described in Beers et al. (1990). The values obtained with this method, together with the best fit values, are listed in Table 1.

The $y$ radial profiles recovered for the cluster subsamples considered are displayed in Fig. 4, together with the joint best fit $y$ profile model obtained, as shown in Table 1.

Finally, the mean masses recovered for the cluster subsamples are listed in Table 4.3.