A COMPARISON AND JOINT ANALYSIS OF SUNYAEV–ZEL'DOVICH EFFECT MEASUREMENTS FROM PLANCK AND BOLOCAM FOR A SET OF 47 MASSIVE GALAXY CLUSTERS

The Sunyaev–Zel'dovich (SZ) effect has emerged as a valuable observational tool for studying galaxy clusters, particularly with the dramatic improvements in instrumentation that have occurred over the past decade. For example, the South Pole Telescope (SPT, Bleem et al. 2015), the Atacama Cosmology Telescope (ACT, Hasselfield et al. 2013), and Planck (Planck Collaboration et al. 2015c) have delivered catalogs with a combined total of more than 1000 SZ-detected clusters. Beyond these large surveys, detailed studies of the gaseous intracluster medium (ICM) have been enabled by an additional set of pointed SZ facilities with broad spectral coverage and/or excellent angular resolution such as the Multiplexed SQUID/TES Array at Ninety GHz (MUSTANG, Mason et al. 2010) the New IRAM KID Arrays (NIKA, Adam et al. 2016), and the Multiwavelength Submillimeter Inductance Camera (Sayers et al. 2016).

As the range of SZ instrumentation has become more diverse, the benefits of joint analyses using multiple data sets have increased. For example, a wide range of studies have used data from two or more SZ receivers in order to measure the spectral shape of the SZ signal (e.g., Kitayama et al. 2004; Zemcov et al. 2010; Mauskopf et al. 2012), mainly for the purpose of constraining the ICM velocity via the kinetic SZ signal, but also to characterize relativistic corrections to the classical SZ spectrum (e.g., Sunyaev & Zeldovich 1980; Nozawa et al. 1998; Chluba et al. 2012). Furthermore, recent analyses have begun to exploit the different angular sensitivities of the SZ facilities in order to obtain a more complete spatial picture of the cluster (e.g., Rodríguez-Gonzálvez et al. 2015; Romero et al. 2015; Young et al. 2015).

In order for these joint analyses to be useful, the various SZ instruments must provide measurements of the SZ signal that are consistent. Historically, this was often not the case, likely due to large systematic errors in the measurements (e.g., see the detailed discussion in Birkinshaw 1999). However, the situation has improved considerably with advances in modern SZ instrumentation, and good agreement has been seen in most recent comparisons (e.g., Mauskopf et al. 2012; Reese et al. 2012; Rodríguez-Gonzálvez et al. 2015; Sayers et al. 2016). Modest inconsistencies do still appear, though they are often the result of assuming different spatial templates when performing the SZ analyses for separate instruments (e.g., Benson et al. 2004; Planck Collaboration et al. 2013b; Perrott et al. 2015). In sum, the systematics that plagued early SZ measurements appear to be largely absent from modern data. This fact, combined with the high degree of complementarity between different SZ facilities, has opened a promising future for detailed cluster studies using multiple SZ data sets.

In this work, we use SZ measurements from Planck and the ground-based receiver Bolocam to study a set of 47 massive clusters. The manuscript is organized as follows. In Section 2, the parametric model used to describe the data is introduced, and in Section 3 the SZ data from Planck and Bolocam are detailed. Section 4 compares the SZ signals measured by Planck and Bolocam, and Section 5 presents the results from joint fits to the two data sets. In Section 6, we use Bolocam SZ measurements to perform a test of the Planck cluster survey completeness, and a summary of the manuscript is given in Section 7.

The thermal SZ effect (Sunyaev & Zel'dovich 1972) describes the Compton scattering of CMB photons with hot electrons in the ICM according to

$$\begin{eqnarray*}{\rm{\Delta }}T(\nu ) & = & f(\nu ,{T}_{{\rm{e}}})y\\ y & = & \displaystyle \int \displaystyle \frac{{k}_{{\rm{B}}}{\sigma }_{{\rm{T}}}}{{m}_{{\rm{e}}}{c}^{2}}{p}_{{\rm{e}}}{dl},\end{eqnarray*}$$

where ΔT(ν) is the observed surface brightness fluctuation in units of CMB temperature at the frequency ν, T_e is the ICM electron temperature, f(ν, T_e) describes the spectral dependence of the SZ signal including relativistic corrections (e.g., Rephaeli 1995; Itoh et al. 1998; Nozawa et al. 1998; Itoh & Nozawa 2004; Chluba et al. 2012), y is the SZ Compton parameter, k_B is Boltzmann's constant, σ_T is the Thompson cross section, m_e is the electron mass, c is the speed of light, p_e is the ICM electron pressure, and dl is along the line of sight. In the absence of relativistic corrections, which are generally small and/or constrained using a spectroscopic X-ray measurement of the value of T_e, the SZ brightness gives a direct measure of the integrated ICM pressure. Therefore, SZ measurements are often used to constrain parametric models of the pressure, such as the generalized Navarro, Frenk, and White (gNFW, Hernquist 1990; Zhao 1996; Navarro et al. 1997) model described in the following section.

2.1. The gNFW Model

Nagai et al. (2007) proposed the use of a gNFW model to describe cluster pressure profiles according to

$$\begin{eqnarray*}&&P(R)=\displaystyle \frac{{P}_{0}}{{(R/{R}_{{\rm{s}}})}^{\gamma }{(1+{(R/{R}_{{\rm{s}}})}^{\alpha })}^{(\beta -\gamma )/\alpha }},\end{eqnarray*}$$

where P(R) is the pressure as a function of radius, P₀ is the normalization factor, R_s is the scale radius, and α, β, and γ control the logarithmic slope of the profile at R ∼ R_s, R ≫ R_s, and R ≪ R_s. Often, the radial coordinates are rescaled to angular coordinates denoted by θ and ${\theta }_{{\rm{s}}}$, and R_s is often recast in terms of a concentration parameter, with

$$\begin{eqnarray*}&&{C}_{500}={R}_{500}/{R}_{{\rm{s}}}={\theta }_{500}/{\theta }_{{\rm{s}}},\end{eqnarray*}$$

and R₅₀₀ denoting the radius where the average enclosed density is 500 times the critical density of the universe. Therefore, for a given value of C₅₀₀, the values of R_s and ${\theta }_{{\rm{s}}}$ are directly related to the cluster mass, M₅₀₀. Furthermore, the normalization is often given in terms of the SZ observable integrated within a specific radius, for example

$$\begin{eqnarray*}&&{Y}_{5{\rm{R}}500}={\int }_{0}^{5\times {\theta }_{500}}y\,\times \,2\pi \theta d\theta .\end{eqnarray*}$$

Nagai et al. (2007) noted that, when P₀ is scaled according to a factor that depends on the cluster's mass and redshift and R_s is recast in terms of C₅₀₀, that a single set of values for α, β, and γ provide an approximately universal description of any cluster's pressure profile. Subsequently, several groups have published different values for these logarithmic slopes based on different samples, data, and analysis techniques (e.g., Arnaud et al. 2010; Plagge et al. 2010; Planck Collaboration et al. 2013a; Sayers et al. 2013; McDonald et al. 2014 and Mantz et al. 2016), and the values given by Arnaud et al. (2010) are the most widely used. The corresponding gNFW shape with C₅₀₀ = 1.18, α = 1.05, β = 5.49, and γ = 0.31 is often referred to as the A10 model.

3.1. Cluster Sample

This study focuses on a set of 47 clusters with publicly available data from Bolocam⁷ and Chandra. Data for 45 of these clusters were published in Czakon et al. (2015), who named that sample the Bolocam X-ray SZ (BoXSZ) sample. Throughout this work, the slightly expanded set of 47 clusters is referred to as the BoXSZ⁺ sample (see Table 1). Based on the Planck MMF3 detection algorithm, 32 BoXSZ⁺ clusters were detected by Planck, with 25 detected at a high enough significance to be included in the Planck cluster cosmology analysis (Planck Collaboration et al. 2015b, 2015c).

Table 1.  Cluster Sample

Cluster	z	R.A.	decl.	M500	${\theta }_{500}$ (CXO)	${\theta }_{500}$ (XMM)	Planck	Bolocam
		hr	deg	1014 M⊙	arcmin	arcmin	SNR	SNR
Abell 2204	0.15	16:32:47	+05:34:32	10.3 ± 1.5	9.2 ± 0.4	8.6 ± 0.2	16.3	22.3
Abell 1689	0.18	13:11:29	−01:20:27	10.5 ± 1.5	7.9 ± 0.4	7.6 ± 0.2	16.7	6.2
Abell 0383	0.19	02:48:03	−03:31:46	4.7 ± 0.8	5.8 ± 0.3	⋯	⋯	9.6
Abell 0209	0.21	01:31:53	−13:36:48	12.6 ± 1.9	7.4 ± 0.4	6.6 ± 0.2	17.1	13.9
Abell 0963	0.21	10:17:03	+39:02:52	6.8 ± 1.0	6.0 ± 0.3	5.8 ± 0.2	8.8	8.3
Abell 1423	0.21	11:57:17	+33:36:39	8.7 ± 2.0	6.5 ± 0.5	5.7 ± 0.2	9.7	5.8
Abell 2261	0.22	17:22:26	+32:07:58	14.4 ± 2.6	7.4 ± 0.4	6.1 ± 0.2	13.5	10.2
Abell 0267	0.23	01:52:42	+01:00:29	6.6 ± 1.1	5.5 ± 0.3	5.0 ± 0.3	5.4	9.6
Abell 2219	0.23	16:40:20	+46:42:29	18.9 ± 2.5	7.9 ± 0.4	6.7 ± 0.1	26.3	11.1
RX J2129.6+0005	0.24	21:29:39	+00:05:17	7.7 ± 1.2	5.6 ± 0.3	4.6 ± 0.3	4.8	8.0
Abell 1835	0.25	14:01:01	+02:52:40	12.3 ± 1.4	6.3 ± 0.3	5.8 ± 0.2	14.4	15.7
Abell 0697	0.28	08:42:57	+36:21:56	17.1 ± 2.9	6.5 ± 0.4	5.6 ± 0.1	18.9	22.6
Abell 0611	0.29	08:00:56	+36:03:25	7.4 ± 1.1	4.7 ± 0.2	4.3 ± 0.2	6.8	10.8
Abell 2744	0.31	00:14:15	−30:23:31	17.6 ± 3.0	6.0 ± 0.5	4.9 ± 0.1	14.1	15.9
MACS J2140.2–2339	0.31	21:40:15	−23:39:40	4.7 ± 0.6	3.9 ± 0.1	⋯	⋯	6.5
Abell S1063	0.35	22:48:44	−44:31:45	22.2 ± 3.4	5.9 ± 0.3	4.8 ± 0.1	20.7	13.6
MACS J1931.8–2635	0.35	19:31:49	−26:34:33	9.9 ± 1.6	4.5 ± 0.2	3.9 ± 0.2	6.1	10.1
MACS J1115.8+0129	0.36	11:15:51	+01:29:54	8.6 ± 1.2	4.2 ± 0.2	3.8 ± 0.2	7.1	10.9
MACS J1532.8+3021	0.36	15:32:53	+30:20:58	9.5 ± 1.7	4.3 ± 0.3	⋯	⋯	8.0
Abell 0370	0.38	02:39:53	−01:34:38	11.7 ± 2.1	4.5 ± 0.3	3.9 ± 0.1	7.6	12.8
MACS J1720.2+3536	0.39	17:20:16	+35:36:22	6.3 ± 1.1	3.6 ± 0.2	3.6 ± 0.2	6.5	10.6
MACS J0429.6–0253	0.40	04:29:36	−02:53:05	5.8 ± 0.8	3.4 ± 0.2	⋯	⋯	8.9
MACS J2211.7–0349	0.40	22:11:45	−03:49:42	18.1 ± 2.5	5.0 ± 0.2	4.1 ± 0.1	11.8	14.7
ZwCl 0024.0+1652	0.40	00:26:35	+17:09:40	4.4 ± 1.6	3.1 ± 0.3	⋯	⋯	3.3
MACS J0416.1–2403	0.42	04:16:08	−24:04:13	9.1 ± 2.0	3.8 ± 0.5	1.5 ± 0.2	4.7	8.5
MACS J0451.9+0006	0.43	04:51:54	+00:06:18	6.3 ± 1.1	3.3 ± 0.2	⋯	⋯	8.1
MACS J0417.5–1154	0.44	04:17:34	−11:54:27	22.1 ± 2.7	4.9 ± 0.2	4.0 ± 0.1	13.3	22.7
MACS J1206.2–0847	0.44	12:06:12	−08:48:05	19.2 ± 3.0	4.7 ± 0.2	3.9 ± 0.1	13.3	21.7
MACS J0329.6–0211	0.45	03:29:41	−02:11:46	7.9 ± 1.3	3.4 ± 0.2	⋯	⋯	12.1
MACS J1347.5–1144	0.45	13:47:30	−11:45:08	21.7 ± 3.0	4.8 ± 0.2	3.8 ± 0.1	11.2	36.6
MACS J1311.0–0311	0.49	13:11:01	−03:10:39	3.9 ± 0.5	2.6 ± 0.1	⋯	⋯	9.6
MACS J2214.9–1400	0.50	22:14:57	−14:00:11	13.2 ± 2.3	3.8 ± 0.2	3.3 ± 0.1	8.3	12.6
MACS J0257.1–2325	0.51	02:57:09	−23:26:03	8.5 ± 1.3	3.2 ± 0.2	2.9 ± 0.1	5.4	10.1
MACS J0911.2+1746	0.51	09:11:10	+17:46:31	9.0 ± 1.2	3.3 ± 0.2	2.9 ± 0.1	5.1	4.8
MACS J0454.1–0300	0.54	04:54:11	−03:00:50	11.5 ± 1.5	3.4 ± 0.2	3.1 ± 0.1	7.1	24.3
MACS J1149.6+2223	0.54	11:49:35	+22:24:04	18.7 ± 3.0	4.0 ± 0.2	3.2 ± 0.1	11.3	17.4
MACS J1423.8+2404	0.54	14:23:47	+24:04:43	6.6 ± 0.9	2.9 ± 0.1	⋯	⋯	9.4
MACS J0018.5+1626	0.55	00:18:33	+16:26:13	16.5 ± 2.5	3.8 ± 0.2	3.1 ± 0.1	8.6	15.7
MACS J0717.5+3745	0.55	07:17:32	+37:45:20	24.9 ± 2.7	4.4 ± 0.2	3.4 ± 0.1	12.8	21.3
MACS J0025.4–1222	0.58	00:25:29	−12:22:44	7.6 ± 0.9	2.8 ± 0.1	⋯	⋯	12.3
MS 2053	0.58	20:56:21	−04:37:48	3.0 ± 0.5	2.1 ± 0.2	⋯	⋯	5.1
MACS J0647.8+7015	0.59	06:47:49	+70:14:55	10.9 ± 1.6	3.2 ± 0.2	2.7 ± 0.1	5.8	14.4
MACS J2129.4–0741	0.59	21:29:25	−07:41:31	10.6 ± 1.4	3.1 ± 0.2	⋯	⋯	15.2
MACS J0744.9+3927	0.70	07:44:52	+39:27:27	12.5 ± 1.6	2.9 ± 0.1	⋯	⋯	13.3
CL J1052.7–1357	0.83	01:52:41	−13:58:06	7.8 ± 3.0	2.1 ± 0.6	⋯	⋯	10.2
MS 1054	0.83	10:56:58	−03:37:33	9.0 ± 1.3	2.3 ± 0.2	⋯	⋯	17.4
CL J1226.9+3332	0.89	12:26:57	+33:32:48	7.8 ± 1.1	2.1 ± 0.1	1.9 ± 0.1	4.9	13.0

Note. From left to right, the columns give: the cluster name, redshift, Chandra R.A. centroid, Chandra decl. centroid, Chandra-derived mass, Chandra-derived ${\theta }_{500}$, XMM-like ${\theta }_{500}$, Planck MMF3 SNR, and Bolocam SNR.

Download table as:  ASCIITypeset image

As detailed in Czakon et al. (2015), the BoXSZ sample lacks a well-defined selection function, and the same is true of the slightly expanded set of BoXSZ⁺ clusters used in this work. To address this issue, Czakon et al. (2015) performed an extensive set of simulations to quantify the effects of this ad hoc selection on the SZ/X-ray scaling relations derived in that work. Although they corrected their results based on these simulations, all of the correction factors were less than the uncertainties due to measurement noise. In addition, Sayers et al. (2013) examined the X-ray morphological characteristics of the BoXSZ sample, and concluded that they closely resemble those of the REXCESS clusters, which were specifically chosen as a good representation of the overall cluster population (Böhringer et al. 2007). However, it should be noted that some of the X-ray morphological classes, such as cool-core and non-cool-core clusters, are not uniformly distributed in redshift within the BoXSZ⁺ sample. In summary, the BoXSZ⁺ cluster selection is not expected to produce any significant biases in the results presented below, with the possible exception of trends versus redshift (see Section 5).

3.2. Planck

The 2015 Planck data release⁸ contains a range of products related to the SZ signal toward clusters, and this analysis utilizes two of those products: (1) the R2.08 cluster catalog created with the MMF3 detection algorithm, which was the baseline catalog for the Planck cluster cosmology analysis (Planck Collaboration et al. 2015b), and (2) the R2.00 all-sky y-maps created based on the MILCA algorithm (Planck Collaboration et al. 2015d), which, as detailed below, show good consistency with the MMF3 measurements for the clusters in the BoXSZ⁺ sample.

The MMF3 catalog provides a two-dimensional probability density function (PDF) for each cluster as a function of ${Y}_{5{\rm{R}}500}$ and ${\theta }_{{\rm{s}}}$ assuming an A10 profile. A constraint on ${Y}_{5{\rm{R}}500}$ can therefore be obtained by marginalizing over ${\theta }_{{\rm{s}}}$, either with or without a prior. As an example of such a prior, the MMF3 catalog provides the values of M₅₀₀ derived from the Planck data, based on a scaling relation calibrated using hydrostatic masses from XMM,⁹ and these values of M₅₀₀ provide a direct constraint on ${\theta }_{{\rm{s}}}$ for an assumed value of C₅₀₀ (see Figure 1).

Zoom In Zoom Out Reset image size

In addition to the MMF3 catalog, the value of ${Y}_{5{\rm{R}}500}$ can also be derived using the all-sky MILCA y-map by fitting an A10 model directly to the map according to the following procedure. First, a prior on the value of ${\theta }_{{\rm{s}}}$ from the XMM-like measurements is used to set the angular size of the model. Then, the three-dimensional model of the cluster is projected to a two-dimensional image with the line-of-sight projection extending to a radial distance of 5 × R₅₀₀. Next, the model is convolved with a 10' full-width half-maxima (FWHM) Gaussian profile to match the point-spread function (PSF) of the MILCA y-map, and binned into square pixels with sides of 333. To compare to this candidate model, the full-sky HEALPix ¹⁰ MILCA y-map data are rebinned into 100' × 100' thumbnails centered on each cluster with identical 333 square pixels (see Figure 1). Next, 1000 random noise maps are generated from the sum of the inhomogeneous noise map and the full-sky homogeneous noise spectrum under the assumption that the noise is Gaussian. From these noise realizations, a variance per pixel is computed, and the inverse of this variance is used as a weighting factor when fitting the A10 model to the data. The fits are performed using the generalized least squares routine MPFITFUN (Markwardt 2009), and the only free parameter in the fits is the overall normalization of the A10 model, ${Y}_{5{\rm{R}}500}$.

The homogeneous noise spectrum of the MILCA y-map is not white, and therefore the per-pixel variance of the random noise maps does not fully describe the data. As a result, the weighting factors used in the fits are in general sub-optimal. This causes the derived parameter uncertainties from the fits to be larger than those from an optimal fit, but it does not produce any bias in the parameter values. However, the parameter uncertainties will, in general, be mis-estimated using this procedure. Consequently, rather than estimating these uncertainties using the per-pixel variance, they are determined using the 1000 noise realizations. Specifically, the best-fit model obtained from the data is added to each of the 1000 noise realizations, all of which are then fit using the same procedure as applied to the real data. For each of these fits, the value of ${\theta }_{{\rm{s}}}$ is varied according to its prior, thus fully including these uncertainties. The spread of values obtained for a given parameter based on these 1000 fits then provides the uncertainty on that parameter.

Based on the above fits, the value of ${Y}_{5{\rm{R}}500}$ obtained from the MMF3 catalog is consistent with the value of ${Y}_{5{\rm{R}}500}$ obtained from the MILCA y-map, with a sample-mean ratio of 1.021 ± 0.023 for the 32 BoXSZ⁺ clusters contained in the MMF3 catalog (see Figure 2). Furthermore, the uncertainty on ${Y}_{5{\rm{R}}500}$ is also consistent between the two, with a sample-mean ratio of 0.967 ± 0.031.¹¹ Therefore, on average, ${Y}_{5{\rm{R}}500}$ values obtained from fits to the MILCA y-maps are equivalent to ${Y}_{5{\rm{R}}500}$ values obtained from the MMF3 catalog.

Zoom In Zoom Out Reset image size

One subtlety in this comparison is the fact that 15 BoXSZ⁺ clusters are excluded because they do not appear in the MMF3 catalog. Due to the complicated degeneracy between ${Y}_{5{\rm{R}}500}$ and ${\theta }_{{\rm{s}}}$, combined with the fact that XMM-like ${\theta }_{{\rm{s}}}$ values are not available for these 15 clusters, it is not possible to obtain accurate upper limits on the value of ${Y}_{5{\rm{R}}500}$ for these MMF3 non-detections based on the publicly available Planck data (e.g., using the Planck MMF3 completeness estimate, which is detailed in Section 6). As a result, there is no way to include these 15 clusters in the comparison between MILCA and MMF3 ${Y}_{5{\rm{R}}500}$ measurements using, for example, survival analysis tools. However, as detailed in Footnote 15, the subsample of 32 BoXSZ⁺ clusters in the MMF3 catalog appears to provide an accurate representation of the ensemble properties of the pressure profiles for the full BoXSZ⁺ sample. Therefore, analyses restricted to the subset of 32 MMF3-detected clusters, such as the comparison described above, should be approximately unbiased.

3.3. Bolocam

The Bolocam fit results from Section 3.3 can be directly compared to the Planck-derived results from Section 3.2, which were based on identical A10 model shapes and XMM-like priors on the value of ${\theta }_{{\rm{s}}}$, along with a nearly identical fitting procedure.¹² The weighted mean ratio between the Bolocam and Planck values of ${Y}_{5{\rm{R}}500}$ obtained from these fits is 1.069 ± 0.030. Given Bolocam's 5% calibration uncertainty, which is common to all of the clusters and therefore acts as a 5% uncertainty on this average ratio, this result indicates consistency (see Figure 2).

Other groups have also compared Planck SZ measurements to ground-based data. For example, Planck Collaboration et al. (2013b) fit A10 models to a set of 11 clusters using XMM priors on ${\theta }_{{\rm{s}}}$ and SZ data from the Arcminute Microkelvin Imager (AMI). They found an average ratio of 0.95 ± 0.05 between the values of ${Y}_{5{\rm{R}}500}$ derived from AMI and Planck, indicating good agreement. A later comparison by Perrott et al. (2015), using AMI observations of 99 clusters, found systematically lower values of ${Y}_{5{\rm{R}}500}$ from AMI relative to Planck. However, the value of ${\theta }_{{\rm{s}}}$ was allowed to float in the fits performed in their analysis, and therefore some or all of the difference in ${Y}_{5{\rm{R}}500}$ values may be a result of using different pressure profile shapes when fitting AMI and Planck. More recently, Rodríguez-Gonzálvez et al. (2015) compared SZ measurements from Planck and the Combined Array for Research in Millimeter-wave Astronomy (CARMA-8) for a set of 19 clusters. Like Perrott et al. (2015), they floated the value of ${\theta }_{{\rm{s}}}$ in their fits; though, unlike Perrott et al. (2015), they obtained consistency, with a CARMA-8/Planck ratio of 1.1 ± 0.4.

Motivated by the good agreement between Planck and Bolocam in measuring the value of ${Y}_{5{\rm{R}}500}$ based on identical A10 profile shapes, the data can be combined to jointly constrain a more general gNFW shape. Specifically, given that Bolocam and Planck are most sensitive to the gNFW shape at large radii, the value of the outer logarithmic slope β is allowed to vary in these fits while the other parameters are fixed to the A10 values. In order to apply these fits to the largest sample possible, namely the full set of 47 BoXSZ⁺ clusters, an external prior on the value of ${\theta }_{{\rm{s}}}$ is required due to the fact that XMM-like priors only exist for 32 BoXSZ⁺ clusters. This ${\theta }_{{\rm{s}}}$ prior is obtained from previously published values of M₅₀₀ derived using data from Chandra, mainly from Sayers et al. (2013) based on the analysis methods detailed in Mantz et al. (2010).¹³ Two BoXSZ⁺ clusters are absent from Sayers et al. (2013), and so the Chandra-derived M₅₀₀ of Abell 1689 is obtained from Mantz et al. (2010) and the Chandra-derived M₅₀₀ of Abell 2744 is obtained from Ehlert et al. (2015).

One subtlety is that the Chandra-derived values of M₅₀₀ are systematically larger than the XMM-like values. In particular, the XMM-like M₅₀₀ values are known to be biased low by ≃30% compared to lensing masses (von der Linden et al. 2014; Planck Collaboration et al. 2015b), while the Chandra M₅₀₀ values described above are ≃10% higher than these same lensing masses (Mantz et al. 2014; Applegate et al. 2016). As a result, the XMM-like values of ${\theta }_{{\rm{s}}}$ are smaller than the Chandra values of ${\theta }_{{\rm{s}}}$, with an average ratio of 1.16 for the 32 BoXSZ⁺ clusters in the MMF3 catalog. Since ${\theta }_{{\rm{s}}}$ sets the angular scale of the gNFW profile, this is equivalent to a change in the value of C₅₀₀. However, since β is allowed to vary in these fits, and β and C₅₀₀ are highly degenerate over the angular scales probed by Planck and Bolocam, this differing choice of ${\theta }_{{\rm{s}}}$ values does not significantly impact the derived profile shape in the radial range where Bolocam and Planck are sensitive. Though, as detailed below, the specific value of β derived from these fits does depend on the choice of ${\theta }_{{\rm{s}}}$ (i.e., of C₅₀₀).

To better understand the results of these jointly constrained gNFW models, linear fits of P₀ and β were performed versus M₅₀₀ and z using LINMIX_ERR (Kelly 2007), which is a Bayesian linear regression routine that allows for measurement errors in both variables and fits for an intrinsic scatter along with the slope and intercept.¹⁴ Uniform priors were assumed for all three fit parameters, and the results of these fits are shown in Figure 3. None of the fits produce strong statistical evidence for a correlation between the measured values, with β versus M₅₀₀ showing the most significant hint of a correlation among the various combinations. These fits find a cluster-to-cluster scatter of ≃30% for the value of P₀ and ≃15% for the value of β. If the linear fits versus mass are evaluated at the median value for the BoXSZ⁺ sample, M₅₀₀ = 9.5 × 10¹⁴ M_☉, then the results are P₀/P₅₀₀ = 9.13 ± 0.68 ± 2.98 and β = 6.13 ± 0.16 ± 0.76 (where the first value represents measurement uncertainty and the second indicates intrinsic cluster-to-cluster scatter). Compared to the A10 model, with P₀/P₅₀₀ = 8.40 and β = 5.49, both of these values are slightly larger and indicate a higher pressure in the cluster center with a steeper fall-off at large radius.

Zoom In Zoom Out Reset image size

However, in comparing these results, it is important to understand how the parameter values are impacted by the choice of mass calibration used to set ${\theta }_{{\rm{s}}}$. Specifically, when fits are performed using the XMM-like values of ${\theta }_{{\rm{s}}}$, which are based on the same mass calibration used by A10 in their analysis and are available for 32 BoXSZ⁺ clusters, P₀/P₅₀₀ increases by an average of 0.97 ± 0.41 and β decreases by an average of 0.46 ± 0.33 (where the second number represents the cluster-to-cluster rms). These shifts are larger than the measurement uncertainties on each parameter, and indicate that the mass calibration can have a noticeable impact. Assuming these results are applicable to the full BoXSZ⁺ sample,¹⁵ they imply that the use of XMM-like ${\theta }_{{\rm{s}}}$ values would have resulted in P₀/P₅₀₀ ≃ 10.10 and β ≃ 5.67. These results indicate that, although the XMM-like ${\theta }_{{\rm{s}}}$ values produce statistically significant shifts in both P₀/P₅₀₀ and β, the qualitative picture compared to the A10 model remains unchanged (i.e., both parameters are still slightly larger compared to that model).

Furthermore, in interpreting these results, it is important to note that, while β provides one metric for understanding the pressure profile at large radius, it does not uniquely describe a single shape due to the strong degeneracies between the gNFW parameters. A more robust metric is the ratio between the integrated SZ signal at R₅₀₀ and at 5 × R₅₀₀,¹⁶ with Arnaud et al. (2010) obtaining Y₅₀₀/${Y}_{5{\rm{R}}500}$ = 0.56. This result can be compared to the value of Y₅₀₀/${Y}_{5{\rm{R}}500}$ = 0.66 ± 0.02 ± 0.10 obtained using the best-fit β = 6.13 derived from our joint Bolocam/Planck fits to the BoXSZ⁺ clusters (see Table 2). Though, as with the value of β, the value of Y₅₀₀/${Y}_{5{\rm{R}}500}$ also depends on the mass calibration used to define ${\theta }_{{\rm{s}}}$. In particular, the best-fit β ≃ 5.67 inferred for an XMM-like mass calibration would result in Y₅₀₀/${Y}_{5{\rm{R}}500}$ ≃ 0.59, much closer to the A10 value.

Table 2.  gNFW Outer Profile Shapes

Analysis	Data Type	β	Y500/${Y}_{5{\rm{R}}500}$
This work	SZ observations	6.13	0.66
This work (XMM-like)	SZ observations	5.67	0.59
Le Brun et al. (2015)	Simulations	4.63	0.63
Ramos-Ceja et al. (2015)	SZ power spectrum	6.35	0.69
Sayers et al. (2013)	SZ observations	3.67	0.28
Planck Collaboration et al. (2013a)	SZ/X-ray observations	4.13	0.48
Battaglia et al. (2012)	Simulations	5.75	0.63
Plagge et al. (2010)	SZ observations	5.5	0.53
Arnaud et al. (2010)	X-ray observations/simulations	5.49	0.56
Nagai et al. (2007)	X-ray observations/simulations	5.0	0.52

Note. Measurements of the outer pressure profile shape in large samples of clusters. The columns show the reference to the analysis, the type of data used in the analysis, the value of β, and the value of Y₅₀₀/${Y}_{5{\rm{R}}500}$. The top row shows the nominal results from this work, based on Chandra-derived masses to set the value of ${\theta }_{{\rm{s}}}$. The second row provides results using XMM-like masses, which have been approximated based on fits to the 32/47 BoXSZ⁺ clusters that have XMM-like masses available (see the text for more details). The difference between these two rows illustrates the potential impact of mass calibration in interpreting scaled pressure profile shapes. In the case of Le Brun et al. (2015), their "median AGN 8.0" fits were used, and were scaled to the median mass of the BoXSZ⁺ sample using their fitting formulae. In the case of Battaglia et al. (2012), their "AGN Feedback Δ = 500" fits were used, and were scaled to the median mass and redshift of the BoXSZ⁺ sample using their fitting formulae. Uncertainties are not available for most analyses, and so they have been omitted.

Download table as:  ASCIITypeset image

As mentioned in Section 2.1, a range of other analyses beyond Arnaud et al. (2010) have also constrained gNFW profiles in large samples of clusters. In particular, several groups have examined these profiles at large radius using either simulations or SZ observations. For example, recent simulations from both Kay et al. (2012) and Battaglia et al. (2012) note a trend of increasing β with redshift, and both Battaglia et al. (2012) and Le Brun et al. (2015) find increasing values of β with increasing mass. Specifically, evaluating the Le Brun et al. (2015) fits at the median mass of the BoXSZ⁺ sample yields β = 4.63 and Y₅₀₀/${Y}_{5{\rm{R}}500}$ = 0.63, the latter indicating an outer profile shape that is consistent with our joint Bolocam/Planck fit. Battaglia et al. (2012) used a parameterization allowing P₀, β, and C₅₀₀ to vary with mass and redshift according to functional forms described by, for example,

$$\begin{eqnarray*}&&\beta \propto {M}^{{b}_{{\rm{M}}}}{(1+z)}^{{b}_{{\rm{z}}}}.\end{eqnarray*}$$

Evaluating their "AGN Feedback Δ = 500" fit at the median mass and redshift of the BoXSZ⁺ sample results in a value of β = 5.75 and Y₅₀₀/${Y}_{5{\rm{R}}500}$ = 0.63, both in relatively good agreement with our joint Bolocam/Planck fits.

Given the good agreement of our results with Battaglia et al. (2012), we also fit an identical functional form to the joint Bolocam/Planck constraints on β, finding b_M = 0.077 ± 0.026 and b_z = −0.06 ± 0.09. These can be compared to the values of b_M = 0.048 and b_z = 0.615 obtained by Battaglia et al. (2012), though some caution is required because the values of C₅₀₀ were not varied in our fits as they were by Battaglia et al. (2012). These results indicate that the trend in mass seen in the Battaglia et al. (2012) simulations is reproduced in our fits, but the trend in redshift is not.

The lack of a redshift trend could be a result of the BoXSZ⁺ sample selection, which is biased toward relaxed cool-core systems at low-z and toward disturbed merging systems at high-z (see Sayers et al. 2013). For example, dynamical relaxation can be defined quantitatively in terms of X-ray image features, specifically the sharpness of the surface brightness peak and the alignment and symmetry of a series of standardized cluster isophotes (Mantz et al. 2015). To be classified as relaxed according to this symmetry, peakiness, and alignment (SPA) criteria, a cluster must exceed the specified thresholds from Mantz et al. (2015) for all three morphological indicators. Based on the SPA criteria, 13 BoXSZ⁺ clusters are relaxed, and these clusters produce a value of β = 6.83 ± 0.37. A similar criteria can be used to define merging clusters within the BoXSZ⁺ sample. Specifically, merging clusters are required to fall below the symmetry and alignment thresholds of the SPA test¹⁷ or contain a radio relic/halo based on the analysis of Feretti et al. (2012) or Cassano et al. (2013). Based on this criteria, 10 BoXSZ⁺ clusters are defined as merging, and these clusters produce a value of β = 5.59 ± 0.61. Therefore, an excess of cool-core clusters at low-z (which have larger values of β on average), and an excess of merging clusters at high-z (which have smaller values of β on average), will artificially introduce a trend of decreasing β with redshift for the BoXSZ⁺ sample.

Other groups have used SZ observations to constrain gNFW profile shapes at large radii. For example, Plagge et al. (2010) fit SZ data from a set of 15 clusters, finding β = 5.5 and Y₅₀₀/${Y}_{5{\rm{R}}500}$ = 0.53. More recently, Planck Collaboration et al. (2013a) used Planck observations of a larger cluster sample to constrain β = 4.13 and Y₅₀₀/${Y}_{5{\rm{R}}500}$ = 0.48 (see Table 2 and Figure 4). Both of these analyses indicate a shallower outer profile than our joint Bolocam/Planck fits; though, some of this difference may be a result of sample selection. Specifically, the Plagge et al. (2010) sample contains clusters with a median redshift of 0.28 and a median mass of M₅₀₀ ∼ 8 × 10¹⁴ M_☉, and the Planck Collaboration et al. (2013a) sample contains clusters with a median redshift of 0.15 and a median mass of M₅₀₀ = 6.3 × 10¹⁴ M_☉. If the parameterization of Battaglia et al. (2012) is used to rescale their gNFW fits to the median mass and redshift of the BoXSZ⁺ sample, then the resulting value of Y₅₀₀/${Y}_{5{\rm{R}}500}$ from both the Plagge et al. (2010) and the Planck Collaboration et al. (2013a) fits is equal to 0.56, closer to our value of 0.66 ± 0.02 ± 0.10. The remaining difference is thought to be related to the mass estimates used in these analyses. In particular, Planck Collaboration et al. (2013a) used XMM-derived masses to set the value of ${\theta }_{{\rm{s}}}$, and, as detailed above, our best-fit Y₅₀₀/${Y}_{5{\rm{R}}500}$ is ≃0.59 when using XMM-like masses.

Zoom In Zoom Out Reset image size

In another recent work, Sayers et al. (2013) obtained, from a joint fit to Bolocam observations of all the clusters in the BoXSZ sample, β = 3.67 and Y₅₀₀/${Y}_{5{\rm{R}}500}$ = 0.28, with an overall profile that noticeably diverges from our joint Bolocam/Planck fit at large radius. This is particularly surprising because the cluster samples are nearly identical, consistent mass estimates were used to set ${\theta }_{{\rm{s}}}$, and the only significant difference is the inclusion of Planck data in our current analysis. Because the Bolocam observations were made from the ground at a single observing frequency, they have reduced sensitivity to large angular scales as a result of both atmospheric fluctuations and primary CMB anisotropies. In contrast, Planck is able to remove CMB anisotropies via its multiple observing channels, and it is not subject to atmospheric fluctuations. Therefore, the Planck data are likely to provide more robust constraints on large angular scales. Though efforts were made in Sayers et al. (2013) to account for the atmospheric and CMB noise, they may be the primary cause of the shallower outer profile found in that work.

Beyond these SZ observations of large samples of individual clusters, Ramos-Ceja et al. (2015) used measurements of the SZ power spectrum on small angular scales from the SPT (Reichardt et al. 2012) to constrain the average pressure profile shape. They found that the A10 model needs to be adjusted to have an outer slope of β = 6.35 ± 0.19 (Y₅₀₀/${Y}_{5{\rm{R}}500}$ = 0.69 ± 0.03) in order to match the SPT measurements. Furthermore, if this value of β is adopted, then their analysis implies little or no evolution in its value as a function of redshift. These results are consistent with our findings.

An accurate characterization of the completeness of the Planck cluster survey is required for cosmological analyses, and the discrepancy between the Planck cluster and CMB power spectrum cosmological results motivates special attention to such a characterization (Planck Collaboration et al. 2015b). The details of how the completeness is estimated are given in Planck Collaboration et al. (2015c) and summarized below. First, a set of clusters based on spherical profiles obtained from simulated clusters (Le Brun et al. 2014; McCarthy et al. 2014) are inserted into both real and simulated Planck maps. The MMF3 algorithm is then applied to these maps, and the probability of detecting a cluster above a given SNR is determined as a function of Y₅₀₀ and ${\theta }_{{\rm{s}}}$ based on a brute force Monte-Carlo, which has been publicly released as part of the MMF3 catalog. Ideally, the accuracy of the completeness function would be tested using a catalog of real clusters with known positions, ${\theta }_{{\rm{s}}}$, and Y₅₀₀. In the absence of such a catalog, Planck Collaboration et al. (2015c) undertook a somewhat less demanding test using the MCXC (Piffaretti et al. 2011) and SPT (Bleem et al. 2015) cluster catalogs, which contain cluster positions and ${\theta }_{{\rm{s}}}$ values, but not Y₅₀₀ values.

The BoXSZ⁺ sample enables a better approximation of the ideal test of the Planck completeness because it has positions, ${\theta }_{{\rm{s}}}$, and Y₅₀₀ estimates for each cluster. Specifically, the positions and ${\theta }_{{\rm{s}}}$ values are obtained from Chandra, the latter rescaled by a factor of 1.16 to account for the average difference between the Chandra and XMM values. This rescaling is required because XMM-derived ${\theta }_{{\rm{s}}}$ values were used to calibrate the Planck completeness. Although it would be better to use the XMM ${\theta }_{{\rm{s}}}$ values for all of the BoXSZ⁺ clusters, they only exist for the clusters detected by the MMF3 algorithm, significantly limiting the value of such a test. In order to obtain Y₅₀₀ estimates from Bolocam, the following procedure is applied. First, the Bolocam value of Y₅₀₀ for each BoXSZ⁺ cluster is generated from A10 model fits to the Bolocam data using the Chandra value of ${\theta }_{{\rm{s}}}$. Next, for the 32 BoXSZ⁺ clusters in the MMF3 catalog, the Planck value of Y₅₀₀ is derived from the MMF3 PDF using the XMM-like value of ${\theta }_{{\rm{s}}}$ in order to mimic the computation of Y₅₀₀ values used in the Planck completeness estimate. The ratio of the Bolocam and Planck Y₅₀₀ values is then fit as a function of M₅₀₀ using LINMIX_ERR (see Figure 5). The results of this linear fit, including the ≃25% intrinsic scatter, are then used to rescale the Bolocam Y₅₀₀ measurements for all of the BoXSZ⁺ clusters. By fitting versus M₅₀₀, this ensures that the mass dependence of the profile shape found in Section 5 is fully included in the conversion from Bolocam to Planck measurements of Y₅₀₀. As part of this rescaling, an additional 5% uncertainty is added to account for the Bolocam flux calibration uncertainty; though, the overall error budget is dominated by the intrinsic scatter in the linear fit.

Zoom In Zoom Out Reset image size

The Chandra and Bolocam values of ${\theta }_{{\rm{s}}}$ and Y₅₀₀, rescaled to mimic the XMM and Planck values as described in the previous paragraph, are then inserted into the Planck SNR = 6 completeness estimate to determine a detection probability for each BoXSZ⁺ cluster (see Figure 6). One subtlety is that the noise in the Planck maps is not uniform over the full sky, and it is therefore necessary to account for this variation when calculating the detection probability for each BoXSZ⁺ cluster. Specifically, this variation is accounted for by comparing the noise rms within the MILCA y-map thumbnail centered on each cluster to the average noise rms within the region of sky satisfying the cuts used for the Planck cluster analysis. In general, the local noise is within 5% of the average, and the most extreme local noise deviation is 12%.

Zoom In Zoom Out Reset image size

The left panel of Figure 6 shows the probability for every BoXSZ⁺ cluster to be detected by the Planck MMF3 algorithm with an SNR > 6. There are no obvious outliers, with Planck detecting all of the clusters with a probability of ∼1 and none of the clusters with a probability of ∼0. To provide a quantitative test, a simulation was performed based on the estimated detection probabilities. For each run of the simulation, a random value was drawn for each BoXSZ⁺ cluster based on the detection probability distribution for that cluster, and the total cumulative number of detections was computed. The simulation was repeated 10,000 times, and the resulting 68% and 95% confidence regions on the cumulative detections are plotted in the right panel of Figure 6. The average number of detections in the simulations is 27.6%, and 16% of the simulations result in fewer than the actual number of clusters detected by Planck, which is 25.

This result provides a more extensive validation of the Planck completeness estimate; though, Planck does detect slightly fewer clusters than expected given the Bolocam Y₅₀₀ measurements. Such a shortfall could partially explain the tension seen between the CMB-derived and cluster-derived cosmological constraints (Planck Collaboration et al. 2015a, 2015b). For example, Planck Collaboration et al. (2015b) quantifies the level of tension in terms of a cluster mass bias, with a value of $(1-b)=0.58$ required to forge agreement. This is smaller than the true mass bias, with $(1-b)\simeq 0.7\mbox{--}0.8$ found from lensing-based mass calibrations (von der Linden et al. 2014; Hoekstra et al. 2015; Planck Collaboration et al. 2015b); though, half of the remaining 10%–20% difference may be due to Eddington bias in those mass measurements (Battaglia et al. 2016). Following this convention, the discrepancy between the predicted and actual number of Planck detections from the BoXSZ⁺ sample can be recast as an effective mass bias. In order for the average number of predicted detections to equal the actual number of 25, the Bolocam Y₅₀₀ measurements would need to be lower by a factor of 0.88 ± 0.11. Based on the Y₅₀₀/M₅₀₀ scaling relation derived in Planck Collaboration et al. (2014), this corresponds to an effective mass bias factor of $(1-b)=0.93\pm 0.06$. This effective bias is multiplicative with the true mass bias, and would bring the Planck cluster results into better agreement with the Planck CMB results.

We fit A10 models to the Planck MILCA y-maps using an XMM-like prior on the value of ${\theta }_{{\rm{s}}}$, obtaining ${Y}_{5{\rm{R}}500}$ values consistent with those determined from the Planck MMF3 algorithm using the same ${\theta }_{{\rm{s}}}$ prior. We also derived ${Y}_{5{\rm{R}}500}$ from ground-based Bolocam observations, finding a Bolocam/Planck ${Y}_{5{\rm{R}}500}$ ratio of 1.069 ± 0.030. This value is consistent with unity given calibration uncertainties and implies that Bolocam and Planck measure consistent SZ signals. Our results are in good agreement with previous comparisons between Planck and the ground-based AMI and CARMA-8 receivers, which yielded similar consistency.

We also performed joint fits to the Bolocam and Planck data, using a gNFW model with the outer logarithmic slope β allowed to vary with the other shape parameters fixed to the A10 values. These fits produce average values of β = 6.13 ± 0.16 ± 0.76 and Y₅₀₀/${Y}_{5{\rm{R}}500}$ = 0.66 ± 0.02 ± 0.10, which are in good agreement with recent simulations for clusters matching the masses and redshifts of the BoXSZ⁺ sample. Compared to simulations, our data are also consistent with the trend of increasing β with increasing cluster mass, but they do not reproduce the relatively strong trend of increasing β with increasing redshift, likely due to selection effects in the BoXSZ⁺ sample. Previous SZ measurements of β and Y₅₀₀/${Y}_{5{\rm{R}}500}$ indicate lower values than our results; though, some or all of this difference may be due to a combination of different median masses and redshifts within those samples, different mass measurements used to set the cluster radial scale, and/or measurement noise.

Using Bolocam measurements of Y₅₀₀ and Chandra measurements of ${\theta }_{{\rm{s}}}$, both rescaled to account for systematic differences relative to Planck measurements of Y₅₀₀ and XMM measurements of ${\theta }_{{\rm{s}}}$, we compute the detection probability for each BoXSZ⁺ cluster using the publicly available Planck completeness estimate. We estimate that Planck should detect an average of 27.6 BoXSZ⁺ clusters above the MMF3 SNR limit for the cosmology sample, a value that is within ≃1σ of the actual number of Planck detections, which is 25. Our results therefore provide a further validation of the Planck completeness estimate. Taking the small discrepancy at face value, however, may suggest that Planck detects fewer clusters than expected. Translated to an effective mass bias, this discrepancy yields $(1-b)=0.93\pm 0.06$. This effective mass bias is multiplicative with the true mass bias determined from lensing measurements, and could partially account for the difference between the Planck cluster-derived and CMB-derived cosmological parameters that has not been explained by the lensing measurements (Planck Collaboration et al. 2015a, 2015b).

We acknowledge the assistance of Kathy Deniston, who provided effective administrative support at Caltech; James Bartlett and Jean-Baptiste Melin, who provided useful discussions; the anonymous referee and AAS statistics editor, who provided useful comments and suggestions. J.S. was supported by a NASA/ADAP award; M.N., G.P., and B.W. were supported by the Caltech Summer Research Connection program. S.R.S. was supported by a NASA Earth and Space Science Fellowship and a generous donation from the Gordon and Betty Moore Foundation.

Facilities: Caltech Submillimeter Observatory - , Planck - European Space Agency's Planck space observatory, Chandra - .