<figure><img src="image_1.png" /></figure>

# Chandra Observations of the Planck Early Sunyaev-Zeldovich Sample: A Reexamination of Masses and Mass Proxies

Felipe Andrade-Santos<sup>1,2</sup>①, Gabriel W. Pratt<sup>3</sup>, Jean-Baptiste Melin<sup>4</sup>, Monique Arnaud<sup>3</sup>, Christine Jones<sup>2</sup>, William R. Forman<sup>2</sup>①, Etienne Pointecouteau<sup>5,6</sup>①, Iacopo Bartalucci<sup>3</sup>①, Alexey Vikhlinin<sup>2</sup>①, Stephen S. Murray<sup>7,14</sup>, Pasquale Mazzotta<sup>2,8</sup>①, Stefano Borgani<sup>9,10</sup>①, Lorenzo Lovisari<sup>2,11</sup>①, Reinout J. van Weeren<sup>12</sup>①, Ralph P. Kraft<sup>2</sup>①, Laurence P. David<sup>2</sup>, and Simona Giacintucci<sup>13</sup>①

## Abstract

Using Chandra observations, we derive the $Y_X$ proxy and associated total mass measurement, $M_{500}^{Y_X}$, for 147 clusters with $z < 0.35$ from the Planck early Sunyaev-Zeldovich catalog, and for 80 clusters with $z < 0.22$ from an X-ray flux-limited sample. We reextract the Planck $Y_{SZ}$ measurements and obtain the corresponding mass proxy, $M_{500}^{SZ}$, from the full Planck mission maps, minimizing Malmquist bias due to observational scatter. The masses reextracted using the more precise X-ray position and characteristic size agree with the published PSZ2 values, but yield a significant reduction in the scatter (by a factor of two) in the $M_{500}^{SZ}-M_{500}^{Y_X}$ relation. The slope is $0.93 \pm 0.03$, and the median ratio, $M_{500}^{SZ}/M_{500}^{Y_X} = 0.91 \pm 0.01$, is within the expectations from known X-ray calibration systematics. $Y_{SZ}/Y_X$ is $0.88 \pm 0.02$, in good agreement with predictions from cluster structure, and implying a low level of clumpiness. In agreement with the findings of the Planck Collaboration, the slope of the $Y_{SZ}-D_A^{-2}Y_X$ flux relation is significantly less than unity ($0.89 \pm 0.01$). Using extensive simulations, we show that this result is not due to selection effects, intrinsic scatter, or covariance between quantities. We demonstrate analytically that changing the $Y_{SZ}-Y_X$ relation from apparent flux to intrinsic properties results in a best-fit slope that is closer to unity and increases the dispersion about the relation. The redistribution resulting from this transformation implies that the best-fit parameters of the $M_{500}^{SZ}-M_{500}^{Y_X}$ relation will be sample-dependent.

*Unified Astronomy Thesaurus concepts: Galaxy clusters (584); Intracluster medium (858); Sunyaev-Zeldovich effect (1654)*

*Supporting material: machine-readable tables*

## 1. Introduction

Galaxy clusters reside in the highest ranges of the halo mass function. In the standard $\Lambda$CDM cosmology, massive halos form by the accretion of smaller subclumps (e.g., Jones et al. 1979; Forman & Jones 1982; Allen et al. 2011; Kravtsov & Borgani 2012). Under the influence of gravity, uncollapsed and collapsed subclumps fall into larger halos and occasionally, objects of comparable mass merge with one another. X-ray observations of substructures in galaxy clusters (see, for instance, Jones & Forman 1984, 1999; Mohr et al. 1995; Buote & Tsai 1996; Jeltema et al. 2005; Böhringer et al. 2010; Laganá et al. 2010; Andrade-Santos et al. 2012, 2013) and measurements of the growth of structures (Vikhlinin et al. 2009b; Mantz et al. 2010; Allen et al. 2011; Benson et al. 2013;

Planck Collaboration XX 2014; Planck Collaboration XXIV 2016) show that these objects are still in the process of formation.

A well-established relation between the cluster mass and its observables (such as X-ray luminosity and gas temperature) is crucial to any work that explores the theoretical relation between the number density of collapsed halos (the mass function) and the underlying cosmological parameters (Vikhlinin et al. 2009b; Planck Collaboration XX 2014; Planck Collaboration XXIV 2016; Pratt et al. 2019). The relations between a cluster's observables and its mass are a direct consequence of gravity, which is the main force driving cluster evolution (Kaiser 1986). Departures from the purely gravitational self-similar expectation are usually attributed to nongravitational processes, such as radiative cooling and active galactic nucleus feedback (e.g., Voit et al. 2005; Pratt et al. 2010).

The Planck satellite surveyed the entire sky across nine microwave frequency bands. The resulting data set allowed the

<footnote n="1"><sup>1</sup> Clay Center Observatory, Dexter Southfield, 20 Newton Street, Brookline, MA 02445, USA; fsantos@cfa.harvard.edu</footnote>

<footnote n="2"><sup>2</sup> Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="3"><sup>3</sup> AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="4"><sup>4</sup> IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="5"><sup>5</sup> CNRS; IRAP; 9 Av. Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France</footnote>

<footnote n="6"><sup>6</sup> Université de Toulouse; UPS-OMP; IRAP; Toulouse, France</footnote>

<footnote n="7"><sup>7</sup> Department of Physics and Astronomy, The Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA</footnote>

<footnote n="8"><sup>8</sup> Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica 1, I-00133 Roma, Italy</footnote>

<footnote n="9"><sup>9</sup> Dipartimento di Fisica dell'Università di Trieste, Sezione di Astronomia, Via Tiepolo 11, I-34131 Trieste, Italy</footnote>

<footnote n="10"><sup>10</sup> INAF, Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131, Trieste, Italy</footnote>

<footnote n="11"><sup>11</sup> INAF—Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Piero Gobetti 93/3, I-40129 Bologna, Italy</footnote>

<footnote n="12"><sup>12</sup> Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="13"><sup>13</sup> Naval Research Laboratory, 4555 Overlook Avenue SW, Code 7213, Washington, DC 20375, USA</footnote>

<footnote>Received 2020 June 30; revised 2021 April 3; accepted 2021 April 6; published 2021 June 15</footnote>

<footnote n="14"><sup>14</sup> Steve Murray passed away on 2015 August 10. Completing the Chandra observations would not have been possible without his invaluable contributions.</footnote>
detection of galaxy clusters through the thermal Sunyaev–Zeldovich (SZ) effect (Sunyaev & Zeldovich 1972). The thermal SZ effect is the shift of the cosmic microwave background (CMB) spectrum toward higher frequencies caused by the inverse Compton scattering of CMB photons by hot electrons in the intracluster medium (ICM). Being proportional to the line-of-sight integral of the ICM pressure, the SZ signal is expected to be closely related to the cluster total mass (e.g., da Silva et al. 2004; Pike et al. 2014), and the amplitude of its signal is independent of the cluster redshift, since CMB photons are all emitted at a constant redshift of $z \sim 1000$. As a consequence, SZ surveys can potentially build unbiased cluster samples, covering higher redshifts than X-ray samples, and are expected to be close to mass-limited.

The second Planck catalog (PSZ2; Planck Collaboration XXVII 2016), derived from the full High Frequency Instrument (HFI) survey of 29 months, detected 1653 cluster candidates. The vast majority (>1200) of these candidates have been confirmed, making the PSZ2 catalog a reference for cluster studies. Among the many quantities provided by the PSZ2 catalog for each cluster, the mass estimate is arguably the most important. Since the Planck masses were determined using XMM-Newton data to calibrate the $Y_{SZ}-Y_X$ relation, comparing Chandra X-ray-derived masses to SZ-derived masses provides an independent and invaluable test for the Planck mass estimates.

Here we use Chandra observations to compute the masses for 147 Planck early SZ (ESZ) clusters at $z < 0.35$. All of the clusters have Planck SZ masses, and we compare these to the X-ray-derived masses from our Chandra analysis. We compare the results with those obtained from observations of an X-ray-selected sample of 80 of the 100 highest-flux X-ray clusters at Galactic latitudes $|b| > 20^\circ$ and $0.025 < z < 0.22$. We also study the $Y_{SZ}-Y_X$ relation. Our aim is to examine the $Y_{SZ}-Y_X$ relation and to compare X-ray and SZ masses for these clusters. We use the SZ-selected clusters as our baseline sample, and the X-ray-selected objects only for comparison. Examination of the impact of our findings on the Planck SZ cosmological constraints is beyond the scope of the present work.

Throughout this paper, we assume a standard $\Lambda$CDM cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.3$, and $\Omega_\Lambda = 0.7$. The variables $M_{500}$ and $R_{500}$ are the total mass within $R_{500}$ and the radius corresponding to a total density contrast $\Delta = 500$ as compared to $\rho_c(z)$, the critical density of the universe at the cluster redshift. The SZ flux is characterized by $Y_{500}^{SZ}$ (in square arcminutes), or simply $Y_{SZ}$. The quantity $D_A^2 Y_{SZ}$ is then the spherically integrated Compton parameter within $R_{500}$ (in square megaparsecs), where $D_A$ is the angular diameter distance of the cluster. All uncertainties are quoted at the $1\sigma$ (68% confidence) level. The natural logarithm is denoted by $\ln$, while log denotes the decimal logarithm.

## 2. SZ-selected and X-Ray-selected Clusters

For our investigation, we make use of the following cluster samples:

1. An SZ-selected sample of clusters derived from the first catalog of 189 objects detected in the first 10 months of the Planck survey, and released in early 2011 (Planck Collaboration VIII 2011). A Chandra X-Ray Visionary Program (PI: Jones) and High Resolution Camera (HRC) Guaranteed Time Observations (PI: Murray) were combined

to form the Chandra–Planck Legacy Program for Massive Clusters of Galaxies. For each of the 163 ESZ Planck clusters at $z < 0.35$, we obtained Chandra exposures sufficient to collect at least 10,000 source counts.<fnref n="15" /> PSZ2 (Planck Collaboration XXVII 2016), derived from the full-mission maps, provided masses for 157 of the original set of 163 ESZ clusters at $z < 0.35$. The remaining six clusters are not in PSZ2, simply because they fall outside the mask used to construct that catalog (i.e., they are close to point sources). We can thus discard those without changing the sample selection function. We use the PSZ2 mass estimates of the resulting systems as they are largely independent of the original detection (discussed in more detail below). Of the resulting 157 clusters, nine are classified as multiple, defined as having more than one object (visually identified in the X-ray images) at a distance less than $10'$ from the main system. These were removed from the following as the presence of more than one system may lead to a boosting of the SZ signal due to confusion in the Planck beam. We also discarded the very nearby cluster PLCKESZ G234.59 + 73.01 (A1367 at $z = 0.02$) as it is too large to allow a reliable background estimate in the Chandra pointing. Figure 1 displays a mosaic of the 147 clusters from the ESZ sample that are used in this work.

2. An X-ray-selected sample of high-flux clusters. Voevodkin & Vikhlinin (2004) compiled a sample of the 52 X-ray-brightest clusters in the local universe by selecting the highest-flux objects detected in the ROSAT All-sky Survey at $|b| > 20^\circ$ and $z > 0.025$—using the HIFLUGCS<fnref n="16" /> catalog as reference. The sample used here is an extension of Voevodkin & Vikhlinin's (2004) approach, where the flux limit in the 0.5–2.0 keV band was lowered to $f_X > 7.5 \times 10^{-12} \text{ erg s}^{-1} \text{ cm}^{-2}$. The full sample contains 100 clusters; the highest-redshift object is at $z = 0.215$. All have Chandra observations. Eighty-three of the clusters from this sample have Planck SZ masses. Three of these were classified as multiple, according to the definition applied to the ESZ sample, so after removing these clusters we obtained a sample of 80 clusters. Of these, 45 are also in the ESZ sample, and 42 are in the HIFLUGCS sample. Figure 2 displays a mosaic of the 80 clusters from the X-ray sample that are used in this work.

The left panel of Figure 3 presents the redshift distribution of both cluster samples. The Planck-detected clusters are clearly more broadly distributed in redshift than the X-ray clusters. This is due to the nature of the selection: for resolved clusters the SZ signal is independent of the redshift of the cluster because it is the CMB that is distorted (the CMB photons originate at the epoch of recombination—from a constant redshift of $z \sim 1000$), while the X-ray-selected clusters constitute a flux-limited sample, which strongly favors X-ray-brighter, lower-redshift clusters. The right panel of Figure 3 presents the mass (derived using the $Y_X$ proxy) distribution of both samples, showing that the X-ray sample has on average lower mass than the Planck ESZ sample.

<footnote n="15">¹⁵ hea-www.cfa.harvard.edu/CHANDRA_PLANCK_CLUSTER</footnote>

<footnote n="16">¹⁶ The Highest X-Ray Flux Galaxy Cluster Sample (Reiprich & Böhringer 2002; Schellenberger & Reiprich 2017).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Mosaic of the 147 clusters from the ESZ sample that are used in this work. Each cluster image was filtered to 0.5–2.0 keV, background-subtracted, exposure map corrected, binned to 4 × 4 pixels, and smoothed to 3 pixels. The image sizes are 3 × R<sub>500</sub> on each side, where R<sub>500</sub> is obtained from M<sup>YX</sup><sub>500</sub>. The images are corrected for surface brightness cosmological dimming with redshift, and are displayed in logarithmic scale. The color table and scale range is the same for all clusters, so that the images would be identical if the clusters obeyed strict self-similarity. The order of clusters follows that from the machine-readable table available online, here in row-major order. Using standard matrix notation, examples are the Bullet cluster in position c<sub>96</sub>, A3411 in position c<sub>81</sub>, and Coma in position c<sub>33</sub>.</figcaption></figure>

# 3. Data Reduction

## 3.1. Chandra Data

### 3.1.1. Event File Processing

Our Chandra data reduction followed the process described in Vikhlinin et al. (2005). We applied the calibration files of CALDB 4.7.2. The data reduction included corrections for the time dependence of the charge transfer inefficiency and gain, and also a check for periods of high background, which were then omitted. Standard blank sky background files and readout artifacts were subtracted. We also detected compact X-ray

sources in the 0.7–2.0 keV and 2.0–7.0 keV bands using CIAO’s `wavdetect` and then masked these sources, as well as extended substructures, before performing spectral and spatial analyses of the cluster emission. For each cluster, we used all available Chandra ACIS (ACIS-I and ACIS-S) observations within 2 Mpc of the cluster center. The median maximum radius of the observations is 1.53 R<sub>500</sub> (with a maximum radius greater than R<sub>500</sub> for 95% of the sample) and the median azimuthal coverage is 0.92 πR<sup>2</sup><sub>500</sub>. We have checked that our results do not depend either on the maximum radius or on the coverage.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Same as Figure 1, but for the 80 clusters from the X-ray sample that are used in this work. Clusters that are also part of the ESZ sample used in this work are presented from the beginning (top left) until (and including) A3667, which is displayed in position c46.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The orange histograms show the X-ray flux-limited sample, while the blue ones correspond to the ESZ cluster sample. Left: Distribution of cluster redshifts in the ESZ and X-ray-selected samples. The ESZ sample extends to higher redshifts than the X-ray flux-limited sample. Right: Distribution of cluster masses in the ESZ and X-ray-selected samples. The X-ray sample extends to lower masses than does the ESZ sample.</figcaption></figure>

### 3.1.2. Emission Measure Profiles

We refer to Vikhlinin et al. (2006) for a detailed description of the procedures we used to compute the emission measure profile for each cluster. We outline here only the main aspects of the method.

We measured the surface brightness profiles in the 0.7–2.0 keV energy band, which maximizes the signal-to-noise ratio (S/N) in Chandra observations for typical cluster gas temperatures. We used the X-ray peak as the cluster center. Readout artifacts and blank-field background (see Section 2.3.3 of Vikhlinin et al. 2006) were subtracted from the X-ray

images, and the results were then exposure-corrected, using exposure maps computed assuming an absorbed optically thin thermal plasma with $kT = 5.0$ keV, abundance = 0.3 solar, and the Galactic column density and including corrections for bad pixels and CCD gaps, which do not take into account spatial variations of the effective area. We subtracted a small uniform component, corresponding to soft X-ray foreground adjustments, if required (determined by fitting a thermal model in the region of the detector field most distant from the cluster center, properly taking into account the expected thermal contribution from the cluster).
<table><caption>**Table 1**

Cluster Masses Using SZ and X-Ray Proxies for the Clusters in the Planck ESZ Sample</caption><thead><tr><th>Cluster</th><th>R.A.</th><th>Decl.</th><th>z</th><th>M<sub>500</sub><sup>Y</sup><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>σ<sub>M<sub>500</sub><sup>Y</sup></sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>M<sub>500</sub><sup>PSX</sup><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>σ<sub>M<sub>500</sub><sup>PSX</sup></sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>M<sub>500</sub><sup>PSZ2</sup><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>σ<sub>M<sub>500</sub><sup>PSZ2</sup></sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>r<sub>max</sub><br>(R<sub>500</sub>)</th></tr></thead><tbody><tr><td>G000.44-41.83</td><td>21:04:18.603</td><td>-41:20:39.36</td><td>0.165</td><td>4.81</td><td>0.47</td><td>4.89</td><td>0.32</td><td>5.30</td><td>0.33</td><td>1.68</td></tr><tr><td>G002.74-56.18</td><td>22:18:39.822</td><td>-38:53:58.47</td><td>0.141</td><td>5.34</td><td>0.51</td><td>4.24</td><td>0.30</td><td>4.41</td><td>0.30</td><td>1.55</td></tr><tr><td>G003.90-59.41</td><td>22:34:27.334</td><td>-37:44:07.88</td><td>0.151</td><td>9.23</td><td>0.89</td><td>7.26</td><td>0.24</td><td>7.19</td><td>0.26</td><td>1.37</td></tr><tr><td>G006.47 + 50.54</td><td>15:10:56.117</td><td>05:44:40.38</td><td>0.077</td><td>7.06</td><td>0.64</td><td>7.19</td><td>0.17</td><td>7.04</td><td>0.20</td><td>1.43</td></tr><tr><td>G006.70-35.54</td><td>20:34:46.912</td><td>-35:49:24.54</td><td>0.089</td><td>5.22</td><td>0.50</td><td>4.40</td><td>0.21</td><td>4.17</td><td>0.23</td><td>1.66</td></tr><tr><td>G006.78 + 30.46</td><td>16:15:46.073</td><td>-06:08:54.61</td><td>0.203</td><td>17.53</td><td>2.01</td><td>16.25</td><td>0.26</td><td>16.12</td><td>0.29</td><td>1.08</td></tr><tr><td>G008.30-64.75</td><td>22:58:48.095</td><td>-34:48:04.62</td><td>0.312</td><td>7.97</td><td>0.75</td><td>7.70</td><td>0.40</td><td>7.75</td><td>0.40</td><td>1.50</td></tr><tr><td>G008.44-56.35</td><td>22:17:45.701</td><td>-35:43:32.55</td><td>0.149</td><td>3.84</td><td>0.39</td><td>4.41</td><td>0.27</td><td>4.79</td><td>0.30</td><td>1.82</td></tr><tr><td>G008.93-81.23</td><td>00:14:19.305</td><td>-30:23:29.33</td><td>0.307</td><td>10.22</td><td>0.93</td><td>9.23</td><td>0.36</td><td>9.84</td><td>0.39</td><td>1.37</td></tr><tr><td>G018.53-25.72</td><td>20:03:30.848</td><td>-23:23:37.54</td><td>0.317</td><td>9.79</td><td>0.93</td><td>8.65</td><td>0.45</td><td>8.99</td><td>0.47</td><td>1.42</td></tr></tbody></table>

**Note.** The columns list the cluster name (ESZ Planck cluster name with the prefix PLCKESZ omitted), R.A., decl., redshift, $Y_X$-derived Chandra mass, $Y_X$-derived Chandra mass uncertainty, PSX Planck mass, PSX Planck mass uncertainty, PSZ2 Planck mass, PSZ2 Planck mass uncertainty, and maximum cluster radius where the emission integral is computed in terms of $R_{500}$ ($R_{500}$ was computed using the $Y_X$-derived Chandra mass).

(This table is available in its entirety in machine-readable form.)

<table><caption>**Table 2**

Cluster Masses Using SZ and X-Ray Proxies for the Clusters in the X-Ray Sample</caption><thead><tr><th>Cluster</th><th>R.A.</th><th>Decl.</th><th>z</th><th>M<sub>500</sub><sup>Y</sup><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>σ<sub>M<sub>500</sub><sup>Y</sup></sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>M<sub>500</sub><sup>PSX</sup><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>σ<sub>M<sub>500</sub><sup>PSX</sup></sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>M<sub>500</sub><sup>PSZ2</sup><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>σ<sub>M<sub>500</sub><sup>PSZ2</sup></sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>r<sub>max</sub><br>(R<sub>500</sub>)</th></tr></thead><tbody><tr><td>G006.47 + 50.54</td><td>15:10:56.117</td><td>05:44:40.38</td><td>0.077</td><td>7.06</td><td>0.64</td><td>7.19</td><td>0.17</td><td>7.04</td><td>0.20</td><td>1.43</td></tr><tr><td>G006.70-35.54</td><td>20:34:46.912</td><td>-35:49:24.54</td><td>0.089</td><td>5.22</td><td>0.50</td><td>4.40</td><td>0.21</td><td>4.17</td><td>0.23</td><td>1.66</td></tr><tr><td>G006.78 + 30.46</td><td>16:15:46.073</td><td>-06:08:54.61</td><td>0.203</td><td>17.53</td><td>2.01</td><td>16.25</td><td>0.26</td><td>16.12</td><td>0.29</td><td>1.08</td></tr><tr><td>G021.09 + 33.25</td><td>16:32:46.854</td><td>05:34:31.61</td><td>0.151</td><td>8.47</td><td>0.77</td><td>8.07</td><td>0.27</td><td>7.79</td><td>0.30</td><td>1.42</td></tr><tr><td>G029.00 + 44.56</td><td>16:02:14.068</td><td>15:58:16.23</td><td>0.035</td><td>2.86</td><td>0.26</td><td>2.64</td><td>0.13</td><td>3.53</td><td>0.06</td><td>1.01</td></tr><tr><td>G033.78 + 77.16</td><td>13:48:52.710</td><td>26:35:31.20</td><td>0.062</td><td>5.67</td><td>0.51</td><td>4.72</td><td>0.14</td><td>4.47</td><td>0.14</td><td>1.53</td></tr><tr><td>G042.82 + 56.61</td><td>15:22:29.473</td><td>27:42:18.76</td><td>0.072</td><td>5.63</td><td>0.51</td><td>4.44</td><td>0.20</td><td>4.08</td><td>0.19</td><td>1.47</td></tr><tr><td>G044.22 + 48.68</td><td>15:58:21.100</td><td>27:13:47.87</td><td>0.089</td><td>12.36</td><td>1.12</td><td>9.44</td><td>0.18</td><td>8.77</td><td>0.20</td><td>1.25</td></tr><tr><td>G046.50-49.43</td><td>22:10:19.489</td><td>-12:10:10.03</td><td>0.085</td><td>5.09</td><td>0.48</td><td>4.51</td><td>0.20</td><td>4.39</td><td>0.20</td><td>1.60</td></tr><tr><td>G048.05 + 57.17</td><td>15:21:12.694</td><td>30:38:00.59</td><td>0.078</td><td>3.92</td><td>0.36</td><td>3.66</td><td>0.23</td><td>3.59</td><td>0.21</td><td>1.77</td></tr></tbody></table>

**Note.** Same as Table 1, but for the X-ray sample.

(This table is available in its entirety in machine-readable form.)

Following these steps, we extracted the surface brightness in narrow concentric annuli ($r_{out}/r_{in} = 1.05$) centered on the X-ray peak and computed the Chandra area-averaged effective area for each annulus (see Vikhlinin et al. 2005 for details on calculating the effective area). To compute the emission measure and temperature profiles, we assumed spherical symmetry. The spherical assumption is expected to introduce only small deviations in the emission measure profile (Piffaretti et al. 2003). Using the modeled deprojected temperature (see Andrade-Santos et al. 2017), effective area, and metallicity as a function of radius, we converted the Chandra count rate in the 0.7–2.0 keV band into the emission integral, $EI = \int n_e n_p dV$, within each cylindrical shell. Tables 1 and 2 list the maximum cluster radius where the emission integral is computed ($r_{max}$) for each cluster. In the ESZ sample, 142 (97%) clusters have $r_{max} > R_{500}$,<fnref n="17" /> and in the X-ray sample, 72 (90%) clusters satisfy this condition (four of the eight clusters that do not satisfy this condition are also in the ESZ sample). We have an average azimuthal coverage within $R_{500}$ of $(0.8 \pm 0.2) \times \pi R_{500}^2$ for the ESZ clusters and $(0.6 \pm 0.2) \times \pi R_{500}^2$ for the X-ray clusters.

We fit the emission measure profile assuming the gas density profile follows that given by Vikhlinin et al. (2006):

$$n_e n_p = n_0^2 \frac{(r/r_c)^{-\alpha}}{(1 + r^2/r_c^2)^{3\beta-\alpha/2}} \frac{1}{(1 + r^\gamma/r_s^\gamma)^{\epsilon/\gamma}} + \frac{n_{02}^2}{(1 + r^2/r_{c2}^2)^{3\beta_2}}, \quad (1)$$

where the parameters $n_0$ and $n_{02}$ determine the normalizations of both additive components. $\alpha$, $\beta$, $\beta_2$, and $\epsilon$ are indices controlling the slope of the curve at characteristic radii given by the parameters $r_c$, $r_{c2}$, and $r_s$; $\gamma$ controls the width of the transition region given by $r_s$. Although the relation given by Equation (1) is based on a classic $\beta$-model (Cavaliere & Fusco-Femiano 1976), it is modified to account for a central power-law-type cusp and a steeper emission measure slope at large radii. In addition, a second $\beta$-model is included, to better characterize the cluster core. For further details on this equation, we refer the reader to Vikhlinin et al. (2006). In the fit to the emissivity profile, all parameters are free to vary. For a typical metallicity of 0.3 $Z_\odot$, the reference values from Anders & Grevesse (1989) yield $n_e/n_p = 1.1995$. Examples of

<footnote n="17">¹⁷ $R_{500}$ defines the radius at an overdensity 500 times the critical density of the universe at the cluster redshift. This quantity was obtained from the X-ray-derived masses.</footnote>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Example X-ray image (top panel) and projected emissivity (lower left), gas density (lower center), and gas mass (lower right) profiles for the cluster PLCKESZ G000.44-41.83. The top panel shows the 0.5–2.0 keV, background-subtracted, exposure map corrected ACIS-I image. Black ellipses correspond to the masked X-ray point sources and the cross corresponds to the cluster center. The lower left panel shows the projected emissivity profile. The solid line shows the emission measure integral of the best fit to the emissivity profile given by Equation (1). The lower center panel shows the gas density profile. The solid line shows the density profile obtained from the emissivity profile given by Equation (1). The lower right panel shows the integrated gas mass profile $M_{gas} (< R)$, with the dashed vertical line indicating $R_{500}$. The dashed lines in the gas density and mass profiles show the 68% confidence range. This is an example of a cluster with moderate data quality, which illustrates the average data quality for the clusters in our samples. The total filtered Chandra exposure is 14 ks.</figcaption></figure>

projected emissivity and gas density profiles are presented in Figure 4.

### 3.1.3. Cluster Mass Estimates

Using the gas mass and temperature, we estimated the total cluster mass from the $M_{500}-Y_{500}^X$ scaling relation of Vikhlinin et al. (2009b),

$$M_{500}^{Y_X} = E^{-2/5}(z)A_{YM} \left[ \frac{Y_X}{3 \times 10^{14} M_{\odot} \text{keV}} \right]^{B_{YM}}, \quad (2)$$

where $Y_X = M_{gas,500} \times kT_X$; $M_{gas,500}$ is computed using the best-fit parameters of Equation (1), and $T_X$ is the measured temperature in the (0.15–1) × $R_{500}$ range. We fit the spectrum in the 0.7–7.0 keV energy range with an absorbed APEC model, fixing NH to the total Galactic value (molecular plus atomic hydrogen contributions—Willingale et al. 2013). The parameters $A_{YM} = (5.77 \pm 0.20) \times 10^{14} h^{1/2} M_{\odot}$, and $B_{YM} = 0.57 \pm 0.03$ (Vikhlinin et al. 2009b). Here, $M_{500}^X$ is the total mass within $R_{500}$, and $E(z) = [\Omega_M(1+z)^3 + (1 - \Omega_M - \Omega_\Lambda)(1+z)^2 + \Omega_\Lambda]^{1/2}$ is the function describing the evolution of the Hubble parameter
with redshift. Using Equation (2), $R_{500}$ is computed by solving

$$M_{500} = 500\rho_c(4\pi/3)R_{500}^3, \quad (3)$$

where $\rho_c$ is the critical density of the universe at the cluster redshift. In practice, Equation (2) is evaluated at a given radius, whose result is compared to the evaluation of Equation (3) at the same radius. This process is repeated in an iterative procedure, until the fractional mass difference is less than 1%. We estimated $1\sigma$ statistical uncertainties in the $Y_{500}^X$-derived masses using Monte Carlo simulations. We also added to the Monte Carlo procedure a $1\sigma$ systematic uncertainty of 9% in the mass determination, as discussed by Vikhlinin et al. (2009b).

## 3.2. Planck Data

### 3.2.1. Published Masses

Derivation of the published Planck masses is extensively described in Planck Collaboration XXVII (2016). The mass was estimated from the relation between $Y_{500}^{SZ}$ and $M_{500}$ (Planck Collaboration XX 2014):

$$E(z)^{-2/3} \left[ \frac{D_A^2 Y_{500}^{SZ}}{10^{-4} \text{ Mpc}^2} \right] = 10^{-0.186} \left[ \frac{h}{0.7} \right]^{-0.211} \times \left[ \frac{M_{500}^{SZ}}{6 \times 10^{14} M_{\odot}} \right]^{1.789}. \quad (4)$$

The SZ size–flux degeneracy was broken through the use of a scaling relation linking the SZ flux $Y_{500}^{SZ}$ to the angular size, $\theta_{500}$, which is then intersected with the two-dimensional posterior probability contours in the ($Y_{500}^{SZ}$, $\theta_{500}$) plane. For the PSZ2 masses, $M_{500}^{SZ,PSZ2}$, the $Y_{500}^{SZ}-\theta_{500}$ prior was obtained from the intrinsic relation between $M_{500}$ and $\theta_{500}$ for a given redshift $z$:

$$\theta_{500} = \theta_* \left[ \frac{h}{0.7} \right]^{-2/3} \left[ \frac{M_{500}^{SZ}}{3 \times 10^{14} M_{\odot}} \right]^{1/3} E^{-2/3}(z) \times \left[ \frac{D_A(z)}{500 \text{ Mpc}} \right]^{-1}, \quad (5)$$

where $\theta_* = 6.997'$.

We recall that Equation (4) was calibrated on a baseline mass proxy relation linking the hydrostatic mass and the X-ray analog of the SZ Compton parameter, $Y_{500}^X$. As described in Planck Collaboration XXVII (2016), the resulting SZ mass $M_{500}^{SZ,PSZ2}$ can be viewed as the hydrostatic mass expected for a cluster consistent with the assumed scaling relation, at a given redshift, and given the measured ($Y_{500}^{SZ}$, $\theta_{500}$) posterior information. Hereafter we will refer to $Y_{500}^{SZ}$ and $Y_{500}^X$ simply as $Y_{SZ}$ and $Y_X$.

### 3.2.2. Reextracted Masses

As discussed in Appendix A of Planck Collaboration XI (2011), a refined measurement of the SZ flux can be obtained by using higher-quality prior information on the position and size of the cluster. Such information is readily provided by our Chandra observations. In the following, we refer to the reextracted masses as $M_{500}^{SZ,PSX}$. The SZ signal was thus reextracted from the six HFI temperature channel maps corresponding to the full Planck

mission survey. We used full resolution maps of HEALPix (Górski et al. 2005)<fnref n="18" /> $N_{side} = 2048$ and assumed that the beams were described by circular Gaussians. We adopted beam FWHM values of 9.66', 7.22', 4.90', 4.92', 4.68', and 4.22' for channel frequencies of 100, 143, 217, 353, 545, and 857 GHz, respectively. Bandpass uncertainties were taken into account in the flux measurement. Uncertainties due to beam corrections and map calibrations are expected to be small, as discussed extensively in Planck Collaboration VIII (2011), Planck Collaboration IX (2011), Planck Collaboration X (2011), and Planck Collaboration XI (2011).

We extracted the SZ signal using multifrequency matched filters (MMFs; Herranz et al. 2002; Melin et al. 2006). These require information on the instrumental beam, the SZ frequency spectrum, and a cluster profile; noise auto- and cross-spectra are estimated directly from the data. We used the MMF in a targeted mode, extracting the SZ flux from a position centered on the X-ray emission peak using the “universal” pressure profile of Arnaud et al. (2010) as a spatial template, in an aperture $\theta_{500}$ corresponding to the $R_{500}$ determined in Section 3.1.3 above. The extraction was achieved by excising a $10^\circ \times 10^\circ$ patch with pixel size 1.72 centered on the X-ray position from the six HFI maps and estimating the SZ flux using the MMF. The profiles were truncated at $5 R_{500}$ to ensure integration of the total SZ signal. The flux and corresponding error were then scaled to $R_{500}$ using the profile assumed for extraction, yielding a refined estimate of $Y_{SZ}$. The resulting refined mass estimate $M_{500}^{SZ}$ PSX was calculated as described above, using the $M_{500}^{SZ}-Y_{SZ}$ relation and the measured ($Y_{SZ}$, $\theta_{500}$) posterior information.

## 3.3. Fitting Procedure

We fit each scaling relation between a pair of observables (X, Y) with a power law of the form

$$Y/Y_0 = A[(X/X_0)]^B \quad (6)$$

assuming no evolution, as expected in the self-similar model. The reason for the latter decision is that the redshift leverage is insufficient to put a strong constraint on the evolution because the maximum redshift of the sample is $z < 0.35$. Figure 5 shows that the data are indeed consistent with no evolution in the current sample. The pivot points $X_0, Y_0$ were $7 \times 10^{14} M_\odot, 6 \times 10^{-5} Mpc^2$, and $3 \times 10^{-3} arcmin^2$ for $M, Y$ [Mpc$^2$], and $Y$ [arcmin$^2$], respectively.

The fit was undertaken using linear regression in the log–log plane, taking the uncertainties in both variables and the intrinsic scatter into account. The X-ray and SZ measurements are independent, so there is no covariance between statistical errors for the scaling laws studied here.

We used the orthogonal Bivariate Correlated Errors and Intrinsic Scatter (BCES; Akritas & Bershady 1996) linear regression method to perform the baseline fits. BCES is widely used by the astronomical community, giving results that may easily be compared with other data sets fitted using the same method. The intrinsic vertical scatter was computed from the quadratic difference between the raw scatter and that expected from the statistical uncertainties, as described in Pratt et al. (2009) and Planck Collaboration XI (2011).

<footnote n="18">¹⁸ http://healpix.jpl.nasa.gov</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Ratio between $Y_{SZ}$ and $Y_X$, plotted as a function of $E(z)$. The white dashed line and gray envelope show the error-weighted mean and associated $1\sigma$ uncertainty, $Y_{SZ}/Y_X = 0.88 \pm 0.02$.</figcaption></figure>

We also fit the data using a Bayesian maximum likelihood estimation approach, with Markov Chain Monte Carlo (MCMC) sampling. Using the `emcee` package developed by Foreman-Mackey et al. (2013), we write the likelihood as described by Robotham & Obreschkow (2015):

$$ \ln L = \sum_{i=1}^{N} \left[ \ln \frac{B^2 + 1}{\sigma_i^2} - \frac{(\ln(Y_i/Y_0) - \ln A - B \ln(X_i/X_0))^2}{\sigma_i^2} \right] \quad (7) $$

$$ \sigma_i^2 = \sigma^2 + B^2 \sigma_{X,i}^2 + \sigma_{Y,i}^2 \quad (8) $$

with the intrinsic scatter, $\sigma^2$, as a free parameter. We used flat priors on the parameters, with ranges [0.6, 1.4], [-1.0, 1.0], and [0.0, 1.0] for B, $\ln A$, and $\sigma$, respectively. For comparison, we also used the `LINMIX` (Kelly 2007) and `LIRA` (Sereno 2016) Bayesian regression packages, with the default priors. The latter, as well as `LINMIX`, allows us to correct for Eddington bias due to the nonuniform distribution of the covariate in the presence of intrinsic scatter (Sereno & Ettori 2015; Sereno 2016). For `LIRA`, we used the option `sigma.XIZ.0='prec.dgamma'` and a single Gaussian component for the corresponding latent variable (we checked that adding an additional component does not change the results). For `LINMIX`, we used the default setting of three Gaussians.

The fit results are detailed in Table 3, where it can be seen that the differences in slope and normalization between methods are always less than the $1\sigma$ statistical error on the parameters. We also note the excellent agreement between `LIRA` and `BCES`, suggesting that the Eddington bias due to the nonuniform distribution in the X variable is not an issue.

Table 3 also gives the mean ratio between the SZ and X-ray quantities. This was computed by maximizing the likelihood of the logarithm of the ratios in log-log space, taking into account the statistical errors and intrinsic scatter. We also computed the

<table><caption>**Table 3**
Best-fit Parameters for Scaling Relations, Obtained with BCES, LINMIX, MCMC, and LIRA Fitting Procedures</caption><thead><tr><th>Relation</th><th>ln A</th><th>B</th><th>σ<sub>ln,i</sub></th></tr></thead><tbody><tr><td>BCES orthogonal</td><td></td><td></td><td></td></tr><tr><td>D<sub>A</sub><sup>2</sup> Y<sub>SZ</sub> - Y<sub>X</sub></td><td>-0.115 ± 0.016</td><td>0.96 ± 0.03</td><td>0.092 ± 0.028</td></tr><tr><td>Y<sub>SZ</sub> - D<sub>A</sub><sup>-2</sup> Y<sub>X</sub></td><td>-0.106 ± 0.016</td><td>0.89 ± 0.01</td><td>...</td></tr><tr><td>M<sub>500</sub><sup>SZ,PSZ</sup> - M<sub>500</sub><sup>JX</sup></td><td>-0.089 ± 0.011</td><td>0.93 ± 0.03</td><td>0.096 ± 0.012</td></tr><tr><td>M<sub>500</sub><sup>SZ,PSX</sup> - M<sub>500</sub><sup>JX</sup></td><td>-0.101 ± 0.009</td><td>0.93 ± 0.02</td><td>0.051 ± 0.016</td></tr><tr><td>LINMIX</td><td></td><td></td><td></td></tr><tr><td>D<sub>A</sub><sup>2</sup> Y<sub>SZ</sub> - Y<sub>X</sub></td><td>-0.118 ± 0.018</td><td>0.95 ± 0.02</td><td>0.099 ± 0.071</td></tr><tr><td>Y<sub>SZ</sub> - D<sub>A</sub><sup>-2</sup> Y<sub>X</sub></td><td>-0.103 ± 0.015</td><td>0.89 ± 0.02</td><td>0.056 ± 0.049</td></tr><tr><td>M<sub>500</sub><sup>SZ,PSZ2</sup> - M<sub>500</sub><sup>JX</sup></td><td>-0.096 ± 0.012</td><td>0.89 ± 0.03</td><td>0.098 ± 0.048</td></tr><tr><td>M<sub>500</sub><sup>SZ,PSX</sup> - M<sub>500</sub><sup>JX</sup></td><td>-0.105 ± 0.011</td><td>0.92 ± 0.02</td><td>0.055 ± 0.039</td></tr><tr><td>MCMC</td><td></td><td></td><td></td></tr><tr><td>D<sub>A</sub><sup>2</sup> Y<sub>SZ</sub> - Y<sub>X</sub></td><td>-0.118 ± 0.017</td><td>0.94 ± 0.02</td><td>0.085 ± 0.028</td></tr><tr><td>Y<sub>SZ</sub> - D<sub>A</sub><sup>-2</sup> Y<sub>X</sub></td><td>-0.103 ± 0.015</td><td>0.88 ± 0.02</td><td>...</td></tr><tr><td>M<sub>500</sub><sup>SZ,PSZ</sup> - M<sub>500</sub><sup>JX</sup></td><td>-0.096 ± 0.013</td><td>0.90 ± 0.03</td><td>0.095 ± 0.012</td></tr><tr><td>M<sub>500</sub><sup>SZ,PSX</sup> - M<sub>500</sub><sup>JX</sup></td><td>-0.107 ± 0.010</td><td>0.91 ± 0.02</td><td>0.050 ± 0.015</td></tr><tr><td>LIRA</td><td></td><td></td><td></td></tr><tr><td>D<sub>A</sub><sup>2</sup> Y<sub>SZ</sub> - Y<sub>X</sub></td><td>-0.118 ± 0.017</td><td>0.96 ± 0.02</td><td>0.053 ± 0.032</td></tr><tr><td>Y<sub>SZ</sub> - D<sub>A</sub><sup>-2</sup> Y<sub>X</sub></td><td>-0.103 ± 0.015</td><td>0.89 ± 0.02</td><td>0.032 ± 0.019</td></tr><tr><td>M<sub>500</sub><sup>SZ,PSZ2</sup> - M<sub>500</sub><sup>JX</sup></td><td>-0.090 ± 0.013</td><td>0.93 ± 0.04</td><td>0.058 ± 0.031</td></tr><tr><td>M<sub>500</sub><sup>SZ,PSX</sup> - M<sub>500</sub><sup>JX</sup></td><td>-0.102 ± 0.010</td><td>0.93 ± 0.02</td><td>0.032 ± 0.016</td></tr><tr><td>Ratio</td><td>Mean value</td><td>σ<sub>ln,i</sub></td><td></td></tr><tr><td>D<sub>A</sub><sup>2</sup> Y<sub>SZ</sub>/Y<sub>X</sub></td><td>0.882 ± 0.015</td><td>0.094<sup>+0.035</sup><sub>-0.027</sub></td><td></td></tr><tr><td>M<sub>500</sub><sup>SZ,PSZ</sup>/M<sub>500</sub><sup>JX</sup></td><td>0.922 ± 0.09</td><td>0.105<sup>+0.012</sup><sub>-0.013</sub></td><td></td></tr><tr><td>M<sub>500</sub><sup>SZ,PSX</sup>/M<sub>500</sub><sup>JX</sup></td><td>0.910 ± 0.09</td><td>0.057<sup>+0.019</sup><sub>-0.015</sub></td><td></td></tr></tbody></table>

**Notes.** Relations are expressed as $Y/Y_0 = A[(X/X_0)]^B$ (Equation (6)). The pivot points $X_0, Y_0$ were $7 \times 10^{14} M_\odot, 6 \times 10^{-5} M\text{pc}^2$, and $3 \times 10^{-3} \text{ arcmin}^2$ for $M$, ($D_A^2 Y_{SZ}, Y_X$) [Mpc$^2$], and ($Y_{SZ}, D_A^{-2} Y_X$) [arcmin$^2$], respectively. For all the relations $\sigma_{\ln,i}$ is the intrinsic scatter in log-log space, except for those of `LIRA`, for which we give the intrinsic scatter with respect to the latent variable, $\sigma_{\ln,i|Y|Z}$, constrained by the method. The intrinsic scatter is not presented for the BCES orthogonal and MCMC fitting procedures of the $Y_{SZ}-D_A^{-2}Y_X$ relation, because for each of those, the raw and statistical scatters were compatible.

ratio by weighting each value of the logarithm ratio by the quadratic sum of the statistical error and the intrinsic scatter. This scatter and the weighted mean were then estimated simultaneously by iteration. Both methods give consistent results.

In Planck Collaboration XXIX (2014), the $Y_{SZ}-Y_X$ relation for the PSZ1 catalog was corrected for Malmquist bias affecting the $Y_{SZ}$ determination (see their Figure 32). The correction was made on a cluster-by-cluster basis, following the approach proposed by Vikhlinin et al. (2009b). The clusters in our study are from the ESZ catalog, detected in the Planck maps corresponding to the first 10 months of the survey. However, here we use the SZ signal from PSZ2, reextracted from the full-mission maps. Although the maps are not fully independent, they correspond to significantly different noise realizations and we thus expect the reextracted SZ signal from the PSZ2 maps to be much less affected by Malmquist bias than the original ESZ detection signal. As a first approximation, in the following we assume that the SZ signal is independent of the detection and neglect residual selection effects from
<figure><img src="image_8.png" /><figcaption>Figure 6. The $Y_{SZ}-Y_X$ relation for the ESZ sample. Left: Correlation between $D_A^2 Y_{SZ}$ and $Y_X$ in units of square megaparsecs, and comparison with previous works. Right: Correlation between $Y_{SZ}$ and $D_A^{-2} Y_X$ in units of square arcminutes. The SZ signal is reextracted in the full-mission maps at the X-ray position and size.</figcaption></figure>

<figure><img src="image_9.png" /></figure>

Malmquist bias. We discuss the impact of this assumption extensively in Section 5.1 and Appendix A.

## 4. Results

### 4.1. The $Y_{SZ}-Y_X$ Relation

The total cylindrical SZ signal is given by the integral of the Compton parameter, $Y_{SZ} = \int y d\Omega$, where $\Omega$ is the solid angle, and the Compton parameter $y$ is given by

$$y = \frac{\sigma_T}{m_e c^2} \int_l n_e k T_e dl, \qquad (9)$$

where $\sigma_T$ is the Thomson cross section, $m_e c^2$ is the electron rest mass energy, $l$ is the distance along the line of sight, and $k$ is the Boltzmann constant. The total SZ signal within $R_{500}$ can also be expressed (in its spherical form) as<fnref n="19" />

$$Y_{SZ} = \frac{\sigma_T}{D_A^2 m_e c^2} \int_V n_e k T_e dV, \qquad (10)$$

where $D_A$ is the angular size distance of the cluster and $V$ is a sphere of radius $R_{500}$.

The X-ray analog of the SZ Compton parameter is defined as $Y_X = M_{gas} \times k T_X$, where $M_{gas}$ is the gas mass given by

$$M_{gas} = \int_V \rho_{gas} dV = \mu_e m_p \int_V n_e dV, \qquad (11)$$

where $\mu_e$ is the mean molecular weight per electron, $m_p$ is the proton mass, and $T_X$ is the X-ray-measured temperature in the (0.15–1.0) × $R_{500}$ radial range.

On the other hand, the mass-weighted temperature is given by

$$T_{mw} = \frac{\int_V n_e T_e dV}{\int_V n_e dV}, \qquad (12)$$

which, when combined with Equations (9) and (11), results in

$$D_A^2 Y_{SZ} = c^* M_{gas} T_{mw}, \qquad (13)$$

where $c^*$ is given by

$$c^* = \frac{\sigma_T}{m_e c^2} \frac{1}{\mu_e m_p} = 1.38 \times \frac{10^{-5} \text{ Mpc}^2}{10^{14} \text{ keV M}_\odot}, \qquad (14)$$

where we used $\mu_e = 1.1706$, which was calculated assuming a fully ionized plasma with metallicity of 0.3 $Z_\odot$, and reference values from Anders & Grevesse (1989).

In practice, when computing $M_{gas}$, we used the fact that the emission measure is proportional to $\langle n_e^2 \rangle$, averaged in radial annuli, to estimate $\langle n_e \rangle^2$. For a uniform spherical gas distribution, $\langle n_e^2 \rangle = \langle n_e \rangle^2$, so one can recover the gas density. However, the reality is different, typically with $\langle n_e^2 \rangle > \langle n_e \rangle^2$. Defining the clumping factor as $Q = \langle n_e^2 \rangle / \langle n_e \rangle^2$, we have $M_{gas}^{obs} = Q M_{gas}^{true}$. By definition $Y_X = M_{gas}^{obs} T_X$, so we can write

$$\frac{D_A^2 Y_{SZ}}{c^* Y_X} = \frac{1}{Q} \frac{T_{mw}}{T_X}. \qquad (15)$$

The $Y_{SZ}-Y_X$ relation is thus a fundamental probe of the consistency of the SZ and X-ray signals. When extracted within the same aperture (here Chandra $R_{500}$), the ratio does not depend on any mass calibration. It depends on clumpiness, internal structure (mass-weighted temperature versus spectroscopic temperature), and the absolute calibration of the X-ray-measured temperature (see, e.g., the discussion in Adam et al. 2017). The $Y_{SZ}-Y_X$ relation is expected to have a slope of 1 (if the clusters are self-similar in shape), and a normalization less than 1 due to the radial dependence of the temperature. In the following, $Y_X$ will refer to the quantity renormalized by the $c^*$ factor (Equation (14)), expressed in units of square megaparsecs.

The left-hand panel of Figure 6 shows the $D_A^2 Y_{SZ}-Y_X$ relation, i.e., the correlation between the cluster physical properties—the intrinsic Compton parameter and $Y_X$. The slope, $B = 0.96 \pm 0.03$, is consistent with unity at slightly more than $1\sigma$, while the normalization is $A = 0.89 \pm 0.09$. The slope of the relation is

<footnote n="19">¹⁹ The Planck Collaboration derived the Compton integral within a sphere, instead of deriving it along the line of sight to infinity (cylindrical integral).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison between X-ray-derived and SZ-derived masses for clusters in the PSZ2 sample. Left: Planck mass published in PSZ2, $M_{500}^{SZ,PSZ2}$, iteratively estimated from the $M_{500}^{SZ}-Y_{SZ}$ relation and the $Y_{SZ}$—size degeneracy curve. The X-ray mass is derived from Chandra data using the $Y_X$ mass proxy. Right: Same as left panel but with SZ signal and masses reextracted at the X-ray position and size, $M_{500}^{SZ,PSX}$. The blue and red lines correspond to the identity and the best-fit power law, respectively. The red regions correspond to the $\pm 1\sigma$ uncertainty and scatter around the best fit. The gray regions correspond to the $\pm 2\sigma$ intrinsic scatter.</figcaption></figure>

consistent with that found by Planck Collaboration XXIX (2014), displayed in Figure 6 as an orange line, within the statistical errors. The relation found by Rozo et al. (2012) using Chandra data, shown in green in Figure 6, has a somewhat shallower slope.

However, the picture changes when the $Y_{SZ}-D_A^{-2}Y_X$ relation is plotted, i.e., when considering the observed SZ flux (in square arcminutes). When plotted this way, the slope is clearly not consistent with unity, $B = 0.89 \pm 0.01$. The effect is similar to, and in agreement with, that observed in Planck Collaboration XXIX (2014), $B = 0.91 \pm 0.02$, although it is at a higher statistical precision (see their Section 7 and Appendix D).

## 4.2. Comparison between SZ-derived and X-Ray-derived Masses

We next investigated the relation between the published Planck masses, $M_{500}^{SZ,PSZ2}$, and those obtained from the Chandra data via the $M_{500}-Y_X$ relation, $M_{500}^X$, as described above in Section 3.1.3. The left-hand panel of Figure 7 compares the published PSZ2 masses to those from Chandra for the 147 ESZ objects. The slope of the relation, $B = 0.93 \pm 0.03$, is less than unity at slightly more than the $2\sigma$ level. This translates into an $M_{500}^{SZ}/M_{500}^X$ ratio that decreases from ~0.99 to ~0.88 between $M_{500} = 2 \times 10^{14} M_\odot$ and $M_{500} = 10^{15} M_\odot$.

The right-hand panel of Figure 7 shows the SZ mass reextracted at the X-ray position and size, $M_{500}^{SZ,PSX}$, compared to the Chandra values. The main effect is a decrease of the intrinsic scatter by about a factor of two (from 11% to 6%). The slope, $B = 0.93 \pm 0.02$, is unchanged, although the uncertainty is smaller as a result of the decreased scatter. As the X-ray mass is slightly higher than the PSZ2 mass, fixing the size to the X-ray value tends to increase the SZ mass on average, while extracting the SZ signal at the X-ray position has the opposite effect (the blind SZ signal is extracted at the position maximizing the S/N). As a result the median SZ to X-ray mass ratio, $M_{500}^{SZ}/M_{500}^Y = 0.91 \pm 0.09$, is unchanged, within the statistical errors.

In the following, we use the $M_{500}^{SZ,PSX}$ value as reference, as it is directly physically related to the $M_{500}^Y$ quantity, and we will refer to it as $M_{500}^{SZ}$.

## 4.3. X-Ray-selected Clusters

Here we consider the effect of adding the X-ray-selected clusters to the relations discussed above. These clusters are discussed only for comparison to the SZ-selected objects.

The top left panel of Figure 8 shows the $M_{500}^{SZ}-M_{500}^Y$ relation for the X-ray-selected clusters only. It can be seen that these objects clearly populate the low-mass part of the plot, with the majority of the systems having $M_{500} < 5 \times 10^{14} M_\odot$. The effect of adding the SZ-selected clusters, as shown in the right-hand panel of Figure 8, is to transform what was a relation with a slope slightly less than unity (i.e., the SZ-selected clusters only) into one with a clear offset in normalization. For the full sample of SZ-selected plus X-ray-selected clusters, the slope $B = 0.97 \pm 0.02$, and the median ratio $M_{500}^{SZ}/M_{500}^Y = 0.90 \pm 0.14$.

To shed light on the above, we can use the corresponding intrinsic $D_A^2 Y_{SZ}-Y_X$ relation, which is shown in the bottom left-hand panel of Figure 8. We observe that as expected, the SZ-selected and X-ray-selected objects are distributed in the $Y_{SZ}-Y_X$ plane in a similar fashion to their distribution in the $M_{500}^{SZ}-M_{500}^X$ plane: the X-ray-selected clusters tend to have higher $Y_X$ values for a given $D_A^2 Y_{SZ}$.

The bottom right-hand panel of Figure 8 shows the $Y_{SZ}-Y_X$ relation expressed in apparent flux ($Y_{SZ}-D_A^{-2}Y_X$). Here we notice that the X-ray-selected clusters follow the distribution of SZ-selected systems, but are typically distributed toward higher values. This is a consequence of the X-ray-selected systems being on average at lower redshift. As a result, when the relation is plotted in intrinsic quantities these systems are redistributed in the $Y_{SZ}-Y_X$ plane via the $D_A^2$ factor. As noted in Planck Collaboration XXIX (2014), the shallower slope in square arcminutes results in an overestimate of the dispersion
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Top left: $M_{500}^{\text{SZ}}-M_{500}^{\text{YX}}$ relation for X-ray-selected clusters only. Top right: $M_{500}^{\text{SZ}}-M_{500}^{\text{YX}}$ relation with SZ-selected and X-ray-selected clusters indicated. Bottom left: $D_A^2 Y_{SZ}-Y_{X}$ relation, with SZ-selected and X-ray-selected clusters indicated. Bottom right: Same as bottom left panel, but for the flux $Y_{SZ}-D_A^{-2}Y_{X}$ relation. The SZ signal is extracted at the X-ray position and size.</figcaption></figure>

about the relation in square megaparsecs as a result of this
redistribution. This is discussed in more detail below.

**5. Discussion**

The slope of the $Y_{SZ}-D_A^{-2}Y_X$ relation, $B = 0.89 \pm 0.01$, is in good agreement with that found by Planck Collaboration XXIX (2014), and points to a very significant departure from unity. This is not the case for the $D_A^2 Y_{SZ}-Y_X$ relation between intrinsic properties, indicating that the $Y_{SZ}/Y_X$ variations are more related to the SZ flux than to the intrinsic Compton parameter. This suggests a measurement origin for the deviation of the slope from unity, rather than a physical variation of $Y_{SZ}/Y_X$ (e.g., with mass). We first discuss possible measurement biases that may explain the slope, starting with a possible residual Malmquist bias.

**5.1. Selection Bias, Scatter, and Covariance**

In the study of scaling relations, the importance of taking
into account measurement biases is well recognized (e.g.,
Angulo et al. 2012). Relations between observed quantities
result from the power-law scaling of each observable with
redshift and mass, the most fundamental cluster properties. The
intrinsic scatter of each quantity around the mean relation may
be correlated—through a common physical origin, for instance.
For a cluster sample constructed from a survey, an important
measurement bias is the so-called Malmquist bias. It is most
critical for the observable used for the detection, $O_{\text{det}}$, when it is
measured from the survey data. At a given mass, clusters with
values that are scattered to higher values of $O_{\text{det}}$, either by noise
or by the intrinsic scatter between the observable and the mass,
will be preferentially detected. The average value of $O_{\text{det}}$ will
then be biased high, particularly when close to the detection
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Simulated $Y_{SZ}$ reextracted from PSZ2 maps, as a function of $D_A^{-2}Y_M$ (left columns) and $D_A^{-2}Y_X$ (other columns), for various assumptions on scatter and covariance. $Y_M$ is the Compton parameter corresponding to the cluster mass for the mean $D_A^{-2}Y_{SZ}-M_{500}$ relation, i.e., the value that would be obtained if there were no scatter in the relation. The scatter in $D_A^{-2}Y_{SZ}$ and $Y_X$ at fixed mass, $\sigma_{\log,YSZ}$ and $\sigma_{\log,Y_X}$, respectively, is indicated in each panel. Middle panels: There is no correlation between $Y_{SZ}$ and $Y_X$ deviations. Bottom right: There is covariance with $r = 0.4$ (Equation (A1)). The different cases are discussed in the text.</figcaption></figure>

threshold, affecting the apparent slope of any scaling relation with $\mathcal{O}_{\text{det}}$. This effect is amplified by the nature of the mass function, the number of clusters increasing with decreasing mass (Allen et al. 2011, see their Figure 5).

One may wonder if the deviation from unity of the slope of the $Y_{SZ}-D_A^{-2}Y_X$ relation is not simply due to residual Malmquist bias. Such a residual bias may be due to correlation between the noise in the ESZ and PSZ2 maps, as well as to the intrinsic scatter of the Compton parameter at a given mass. The intrinsic scatter of the $Y_X$ proxy with the underlying mass should also be considered, as should its probable correlation with the scatter in $Y_{SZ}$.

Appendix A describes extensive simulations that were designed to investigate these effects on the observed $Y_{SZ}-D_A^{-2}Y_X$ relation. These were based on the production of mock ESZ catalogs through the injection of simulated clusters, following the Tinker mass function, into the Planck 10 month maps. This was followed by the application of the MMF algorithm to obtain SZ detections above an S/N = 6, the threshold of the ESZ catalog. The detected clusters were then reinjected into the full-mission maps and reanalyzed by the MMF algorithm to extract their PSZ2 flux values. Ten mock catalogs were produced, with a total of 1657 detections. To investigate selection, scatter, and covariance effects, we assumed a Gaussian lognormal correlated distribution for $D_A^{-2}Y_{SZ}$ and $Y_X$ at fixed mass $M_{500}$ with a covariance matrix. The

reference quantity at a given mass, exactly following the scaling relation, was denoted by $Y_M$. The $Y_{SZ}-D_A^{-2}Y_X$ relation of the mock data was then analyzed exactly as that for the real data. Introducing each effect successively (see Figure 9), we obtained the following findings:

1. The residual Malmquist bias, due to the fact that the noise realizations of the ESZ and PSZ2 maps are not strictly fully independent, is negligible. The $Y_{SZ}$ flux extracted from the PSZ2 map is unbiased as compared to the reference $D_A^{-2}Y_M$ value, in the absence of scatter. This is also confirmed by the absence of correlation, in the real data, between $Y_{SZ}/Y_X$ and S/N of the ESZ detection at a given mass (Figure 10). The correlation coefficients for the full sample and subsample are given in the plots of Figure 10 and in Table 4.

2. There is a small residual Malmquist bias due to the scatter in the $D_A^{-2}Y_{SZ}-M_{500}$ relation. This is expected, as clusters detected above the ESZ survey threshold, because they are scattered up due to the intrinsic scatter, will still have a higher $Y_{SZ}$ signal than the average in the PSZ2 maps. The importance of this effect in the $Y_{SZ}-D_A^{-2}Y_M$ relation depends on the relative magnitude of the ESZ noise and the intrinsic scatter, i.e., on which effect drives the cluster detectability. For reasonable values of scatter derived from numerical simulations of cluster formation, the effect on the slope is modest, $\Delta B \sim -0.02$.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Verification of the importance of the Malmquist bias. Left panel: Correlation between $M_{500}^{SZ}$ and the S/N of the ESZ detection. Red points correspond to the subsample used to check for the Malmquist bias effect, effectively removing the intrinsic mass dependence from the S/N of the detection (see text). Center and right panels: Correlation between the SZ and mass ratios as a function of S/N. The correlation coefficients for the full sample and subsample are given in the plots and in Table 4.</figcaption></figure>

<table><caption>**Table 4**
Spearman Rank Coefficient and Significance of Malmquist Bias Tests</caption><thead><tr><th>Relation</th><th>r</th><th>p</th><th>r subsample</th><th>p subsample</th></tr></thead><tbody><tr><td>M<sub>YX</sub>-S/N</td><td>0.542</td><td>9.2 × 10<sup>-13</sup></td><td>-0.003</td><td>0.99</td></tr><tr><td>Y<sub>SZ</sub>/Y<sub>X</sub>-S/N</td><td>-0.345</td><td>1.6 × 10<sup>-5</sup></td><td>-0.015</td><td>0.92</td></tr><tr><td>M<sub>500</sub><sup>SZ</sup>/M<sub>500</sub><sup>X</sup>-S/N</td><td>-0.372</td><td>3.0 × 10<sup>-6</sup></td><td>-0.025</td><td>0.87</td></tr></tbody></table>

**Notes.** Columns list the rank coefficient and significance of Malmquist bias tests for the full sample and subsample. The null hypothesis *p*-value is the probability that the observed coefficient is obtained by chance if the two parameters are completely independent.

3. Intrinsic scatter in the $Y_X-M_{500}$ relation decreases the slope of the $Y_{SZ}-D_A^{-2}Y_X$ relation, as compared to the $Y_{SZ}-D_A^{-2}Y_M$ relation, by $\Delta B \sim -0.02$. We do not have a clear explanation for this effect, which cannot simply be due to a residual Eddington-like bias. However the likely correlation between the $D_A^2 Y_{SZ}$ and $Y_X$ scatter redresses the slope back toward unity ($B = 0.973 \pm 0.005$). This is a well-known effect: selection bias due to intrinsic scatter has little impact on the relation between observables when clusters scatter almost along the relation (Angulo et al. 2012).

In summary, residual Malmquist bias due to Planck noise (statistical scatter) is negligible, and selection effects due to typical values of intrinsic scatter as predicted by numerical simulations cannot explain the observed behavior of the $Y_{SZ}-D_A^{-2}Y_X$ relation.

An independent mass estimate (such as a lensing mass) is required to make further progress in this issue. This would enable one to constrain simultaneously the $Y_X-M_{500}$ and $D_A^2 Y_{SZ}-M_{500}$ scaling relations and their intrinsic scatter, taking into account the Planck selection function. Such an approach, based on a full modeling of the maximum likelihood of the observed data ($Y_{SZ}$, $Y_X$, and mass for each cluster), is now frequently applied to survey data (e.g., Stanek et al. 2006; Vikhlinin et al. 2009b; Mantz et al. 2010; Giles et al. 2016; Dietrich et al. 2019). Application of such a method to the current data set would be more complicated, however, as we would need to

combine the detection in the ESZ maps with the measurement of the $Y_{SZ}$ flux in the PSZ2 maps. Such a study is outside the scope of the paper.

## 5.2. Other Measurement Systematic Effects

Following the tests performed in Section 5.1 and in Appendix A, indicating that selection effects cannot explain the slope of the $Y_{SZ}-D_A^{-2}Y_X$ relation, we next considered the possibility that its departure from unity is due to systematic issues with the $Y_{SZ}$ or $Y_X$ estimates. We therefore considered a range of possible explanations, including the following:

1. Relativistic effects: The results presented in this work neglect relativistic effects, which might have an impact on the measured fluxes. In principle the relativistic effect is stronger in hotter, more massive clusters, so it should be most apparent in the slope of the correlation between intrinsic properties. This is clearly not the case (Figure 6). However, to further test this hypothesis, we reextracted Planck fluxes, assuming relativistic corrections to the SZ spectrum (Itoh et al. 1998) based on the observed Chandra temperatures. Refitting the data obtained with these relativistic corrections, we found the slope of the $Y_{SZ}-D_A^{-2}Y_X$ relation remains identical ($0.89 \pm 0.01$). The slope of the $D_A^2 Y_{SZ}-Y_X$ relation increases only slightly, to $0.99 \pm 0.03$. We conclude that relativistic corrections cannot explain the deviation of the slope from unity for the observed $Y_{SZ}-D_A^{-2}Y_X$ relation.

2. Non-Gaussian beams: The approximation of the Planck beam with a Gaussian may cause a loss of flux due to neglect of the side lobes. This effect would be more important for the most extended objects (i.e., nearby or hot clusters). To test the magnitude of this effect, we reextracted the $Y_{SZ}$ values using the Planck effective beams available through the Planck Legacy Archive.<fnref n="20" /> We found a median ratio of $Y_{SZ}^{\text{Gauss}}/Y_{SZ}^{\text{eff.}} = 1.00 \pm 0.02$, with essentially zero dependence on angular extent.

3. Variation in profiles: As detailed in Section 3.2, reextraction of the Planck flux was undertaken assuming the “universal” pressure profile of Arnaud et al. (2010) as

<footnote n="20">²⁰ https://pla.esac.esa.int/</footnote>
a fixed spatial template. It is well known that individual profiles may deviate significantly from this. Furthermore, nearer clusters are better resolved in the Planck maps than more distant objects, so there may be a distance dependence in this variation. However, comparison of $Y_{SZ}$ fluxes measured using a fixed “universal” template and those based on the individual pressure profile yielded a best-fit slope and normalization entirely consistent with unity (see Figure B.1 of Planck Collaboration Int. V 2013).

4. X-ray temperature calibration issues: The observed good agreement in slope seen in Figure 6 between the present study and the previous work by Planck Collaboration XXIX (2014), which used XMM-Newton data, suggests, empirically, that this cannot explain the observed slope. Furthermore, as for the relativistic effect, any temperature-dependent systematics should be more visible in the correlation in square megaparsecs than in the correlation in square arcminutes, which is clearly not the case. The effect of XMM-Newton versus Chandra calibration on the scaling relations is further discussed in Section 5.4.

In summary, none of the above possibilities explains why the slope of the $Y_{SZ}-Y_X$ relation is significantly less than unity when expressed in units of apparent flux (square arcminutes).

### 5.3. Link between the $M_{500}^{SZ}-M_{500}^{YX}$ and $Y_{SZ}-Y_X$ Relations

The slope of the intrinsic $D_A^2 Y_{SZ}-Y_X$ relation is significantly different from that of the $Y_{SZ}-D_A^{-2}Y_X$ relation, and is consistent with unity. The difference between the two relations is a conversion from flux to intrinsic properties, i.e., a linear rescaling of each quantity by the angular distance $D_A^2$ of each cluster. In Appendix C, we demonstrate analytically that this transformation results in a $D_A^2 Y_{SZ}-Y_X$ relation slope that is closer to unity, as observed ($B = 0.96 \pm 0.03$), and increases the dispersion about the relation. This effect is actually caused by the reordering of the data from the $Y_{SZ}-D_A^{-2}Y_X$ plane to the $D_A^2 Y_{SZ}-Y_X$ plane, as already noted in Section 4.3 for X-ray clusters.

The properties of the $M_{500}^{SZ}-M_{500}^{YX}$ relation, derived from the mass proxies $M_{500}^{SZ}$ and $M_{500}^{YX}$, originate in the properties of the intrinsic $D_A^2 Y_{SZ}-Y_X$ relation. One can translate between the two by making use of the $M_{500}-Y_X$ and $M_{500}-D_A^2 Y_{SZ}$ relations. We can express the relations as follows:

$$ 
\begin{aligned}
D_A^2 Y_{SZ} &\propto (Y_X)^\alpha \\
M_{500}^{SZ} &\propto (M_{500}^{YX})^\beta \\
M_{500}^{YX} &\propto (Y_X)^\gamma \\
M_{500}^{SZ} &\propto (D_A^2 Y_{SZ})^\delta,
\end{aligned} 
\quad (16) $$

where $\alpha, \beta, \gamma$, and $\delta$ denote the slope of the relations above (only here). The above leads to $\beta = \alpha \times \delta/\gamma$. In this work we obtained $\alpha = 0.96 \pm 0.03$. From Equation (4) we have $\delta = 1/1.789 = 0.56$, and from Vikhlinin et al. (2009b) we have $\gamma = 0.57$, which leads to $\beta = (0.96 \pm 0.03) \times (0.56)/0.57 = 0.94 \pm 0.03$, which is basically what one gets in Figure 7 from the direct fit. In other words, the deviation of the slope from unity is simply slightly amplified when considering the mass proxies, rather than the measured $Y$ parameters, because the $M_{500}-D_A^2 Y_{SZ}$ and $M_{500}-Y_X$ relations do not have exactly the same slopes (although they are consistent).

Furthermore, the logarithmic intrinsic scatter of $M_{500}^{SZ}-M_{500}^{YX}$ is $\sigma_{\ln,i} = 0.051$, about $\gamma$ times the intrinsic scatter of the $D_A^2 Y_{SZ}-Y_X$ relation, $\sigma_{\ln,i} = 0.092$, as expected.

### 5.4. Effect of X-Ray Temperature Calibration and of a Different $M_{500}-Y_X$ Calibration

It is well known that there is a systematic difference between the X-ray temperatures derived from Chandra and XMM-Newton, with the former yielding temperatures ~7% higher than XMM-Newton in the mass range considered here (Schellenberger et al. 2015). Consistently, there is also a slight difference of ~4% between the normalization of the $M_{500}-Y_X$ relation obtained by Vikhlinin et al. (2009b), which is used here, and that of the relation obtained by Arnaud et al. (2010). (The slope difference between the two relations of $0.01 \pm 0.02$ is completely negligible.) Here we estimate analytically the effect of using XMM-Newton data and the Arnaud et al. (2010) relation on the $Y_{SZ}-Y_X$ and $M_{500}^{SZ}-M_{500}^{YX}$ relations.

#### 5.4.1. Effect on Measured Quantities

We consider first the effect on $Y_X$ and the corresponding mass $M_{500}^{YX}$. We write the relation as $M_{500}^{YX} \propto AY_X^B$, with $Y_X = M_g(r < R_{500})T_X$. The quantities $M_{500}^{YX}$ and $Y_X$ are estimated iteratively, and $Y_X$ does not simply scale linearly with $T_X$ as $M_g$ depends sensitively on the aperture. For simplicity, we assume that the gas density is given by a $\beta$-model with $\beta = 2/3$, which leads to a gas mass $M_g(r < R_{500}) \propto R_{500}$. On the other hand, we can consider that the gas masses derived from XMM-Newton and Chandra data are the same (Bartalucci et al. 2017) and the core-excised gas temperature $T_X$ is not changed by a (small) change in aperture. Noting that $M_{500}^{YX} \propto R_{500}^3$ by definition, we then obtain

$$ 
\begin{aligned}
Y_X &\propto A^{1/(3-B)} T_X^{3/(3-B)} = A^{0.41} T_X^{1.23} \\
M_{500}^{YX} &\propto A^{3/(3-B)} T_X^{3B/(3-B)} = A^{1.23} T_X^{0.70} \\
R_{500} &\propto A^{1/(3-B)} T_X^{B/(3-B)} = A^{0.41} T_X^{0.23}
\end{aligned} 
\quad (17) $$

for $B = 0.57$ (Equation (2)).

The SZ quantities, $Y_{SZ}$ and the corresponding $M_{500}^{SZ}$ reextracted within the X-ray aperture, will also be affected by the $R_{500}$ aperture change. From our analysis of Planck data for the present sample, the variation of $Y_{SZ}$ with aperture $\theta_{500}$ around the nominal $\theta_{500}$ value has a logarithmic slope of 0.8, so that

$$ Y_{SZ} \propto A^{0.8/(3-B)} T_X^{0.8 \times B/(3-B)}. \quad (18) $$

The corresponding change in SZ mass, $M_{500}^{SZ}$, can be obtained from the $M_{500}-D_A^2 Y_{SZ}$ relation (Equation (5)). We can then insert the numerical values for the Chandra and XMM-Newton temperature difference and for the normalization of the $M_{500}-Y_X$ relation: $T_{X,CXO} = 1.07 T_{X,XMM}$; $A_{XMM} = 0.96 A_{CXO}$ with $B = 0.57$. We then obtain $Y_{X,XMM} \sim 0.90 Y_{X,CXO}$, $M_{500,XMM}^{YX} = 0.88 M_{500,CXO}^{YX}$, $Y_{SZ,XMM} \sim 0.97 Y_{SZ,CXO}$, and $M_{500,XMM}^{SZ} \sim 0.99 M_{500,CXO}^{SZ}$.

#### 5.4.2. Effect on the Scaling Relations

We can now estimate how these changes in cluster parameters would affect the $D_A^2 Y_{SZ}-Y_X$ and $M_{500}^{SZ}-M_{500}^{YX}$ relations. For constant multiplying factors that change $Y_X, Y_{SZ}$, and the corresponding masses, the slopes would not change; however, the normalizations
would. We consider the following power laws: $D_A^2 Y_{SZ} = A(Y_X)^\alpha$ and $M_{500}^{SZ} = B(M_{500}^{Y_X})^\beta$, with $\alpha = 0.96$ and $\beta = 0.93$ corresponding to the best-fitting values (Table 3). The normalization will change following the ratio of the XMM-Newton and Chandra quantities, given above, as $A_{XMM} = A_{CXO} \times [0.97/0.90]^\alpha$, and $B_{XMM} = B_{CXO}[0.99/0.88]^\beta$. Had we used XMM-Newton data and the Arnaud et al. (2010) $M_{500}-Y_X$ relation, we estimate that the slopes for the $Y_{SZ}-Y_X$ and $M_{500}^{SZ}-M_{500}^{Y_X}$ relations would be the same, and the normalizations would be higher, by ~7% and 12%, respectively.

Finally, we compare our results with those from Planck Collaboration XXIX (2014). For comparison purposes, we express the $D_A^2 Y_{SZ}-Y_X$ relation using the pivot point of Planck Collaboration XXIX (2014), $Y_p = 10^{-4}$ Mpc$^2$. In this work, we obtain $Y_{SZ}/Y_p = 0.87 \times (Y_X^{CXO}/Y_p)^{0.96}$. Using the estimated changes in slope and normalization for XMM-Newton data derived above, $\alpha_{XMM} = \alpha_{CXO}$ and $A_{XMM} = 1.07 A_{CXO}$, we obtain $Y_{SZ}/Y_p = 0.93 \times (Y_X^{XMM}/Y_p)^{0.96}$, which is in perfect agreement with the result from Planck Collaboration XXIX (2014), $Y_{SZ}/Y_p = (0.94 \pm 0.02) \times (Y_X/Y_p)^{0.98 \pm 0.03}$. The Planck SZ masses were calibrated from the $Y_{SZ}-Y_X$ and $M_{500}-Y_X$ relations derived from XMM-Newton and the corresponding normalization of the $M_{500}^{SZ}-M_{500}^{Y_X}$ relation is one by construction. The fact that we find a mean mass ratio (0.91 ± 0.09) 10% lower than unity is also consistent with the calibration differences.

In conclusion the impact of XMM-Newton versus Chandra calibration is small, but larger than the statistical errors on mean quantities. There is good agreement in slope between the present study and the previous study by Planck Collaboration XXIX (2014) based on XMM-Newton data. Differences in mean Y and mass ratios are consistent between XMM-Newton and Chandra calibration.

## 5.5. Comparison to Previous Results

Schellenberger & Reiprich (2017) compared the Planck masses with masses derived assuming hydrostatic equilibrium (HE) for 50 clusters that the HIFLUGCS catalog and PSZ2 have in common. They derived a larger discrepancy between X-ray and SZ masses, with a very significant mass dependence (slope of $\alpha = 0.76 \pm 0.08$). They suggested that the differences at low masses can be due to the Malmquist bias effect, in conjunction with an underestimate of the Planck masses at the high-mass end. Our study does not confirm this trend. We emphasize that we have shown that our study is free from Malmquist bias. On the other hand, the selection function of the subsample of HIFLUGCS clusters discussed in Schellenberger & Reiprich (2017) is a complex combination of X-ray and SZ selections. Moreover, their masses were estimated from the HE equation and all the most massive clusters in their sample are highly unrelaxed objects. While generally we expect the HE mass to be lower than the true mass, the HE mass may actually overestimate the true mass when it is obtained from extrapolation of the Navarro–Frenk–White profile fitted to data in the central region (see Figure 5 of Rasia et al. 2006). Lovisari et al. (2020) also compared hydrostatic masses to the Planck masses. They compared the hydrostatic masses of 117 clusters from the ESZ sample with their PSZ2-derived masses, finding a small offset (i.e., 4%) when the enclosed hydrostatic masses were computed at $R_{500}$ as determined by Planck.

# 6. Conclusions

Using Chandra observations, we derived the $Y_X$ proxy and associated mass for 147 clusters with $z < 0.35$ from the Planck ESZ catalog and for 80 clusters with $z < 0.22$ from an X-ray flux-limited sample. The Chandra ESZ follow-up was completed within the Planck cosmological mask region. We reextracted the Planck $Y_{SZ}$ measurements centered on a more precise X-ray position and within the X-ray characteristic size ($R_{500}$) from the full Planck mission maps. This reextraction of the SZ signal for clusters originally detected in the 10 month maps minimized the Malmquist bias, with the residual bias due to statistical scatter being negligible.

We discuss the results in terms of two $Y_{SZ}$ measurements: the apparent flux $Y_{SZ}$ in units of square arcminutes, and the intrinsic Compton parameter $D_A^2 Y_{SZ}$ in units of square megaparsecs. Our conclusions are as follows:

1. The $D_A^2 Y_{SZ}-Y_X$ relation is consistent in slope with that derived by the Planck Collaboration from XMM-Newton data, with $B = 0.96 \pm 0.03$. There is a slight offset in normalization of $6.5\% \pm 0.1\%$ as compared to the Planck results, which, as we have shown, is consistent with known calibration systematic uncertainties between XMM-Newton and Chandra. The resulting $Y_{SZ}/Y_X$ is $0.88 \pm 0.02$, in good agreement with X-ray expectations given radially decreasing temperature profiles, and with previous determinations. It also suggests that there is a low level of gas clumping within $R_{500}$. This result is independent of any mass calibration or scaling relations.

2. We compared the Chandra X-ray masses, derived from the $Y_X$ proxy, to the Planck SZ masses, derived from the $D_A^2 Y_{SZ}$ proxy. The Planck masses are about 10% smaller. The median ratio of $M_{500}^{SZ}/M_{500}^{Y_X} = 0.91 \pm 0.01$ is consistent with the aforementioned calibration differences. The slope of the relation, $0.93 \pm 0.03$, differs from unity at slightly more than the $2\sigma$ level. The use of the reextracted masses does not change the $M_{500}^{SZ}-M_{500}^{Y_X}$ relation, but results in a significant reduction in the scatter (by a factor of two).

3. The slope of the $Y_{SZ}-D_A^{-2}Y_X$ flux relation, at $B = 0.89 \pm 0.01$, is significantly less than unity. We have shown that this effect is not due to measurement issues. We performed extensive simulations, involving injection of simulated clusters into Planck maps, and subsequent detection and reextraction using the Planck SZ detection algorithm MMF3. Using these simulations, we showed that this result is not due to selection effects, intrinsic scatter, or covariance between quantities. The X-ray-selected sample follows the general trend exhibited by the SZ-selected sample, and extends it down to lower SZ fluxes.

4. We have shown analytically that changing the $Y_{SZ}-Y_X$ relation from apparent flux (in units of square arcminutes) to intrinsic properties (in units of square megaparsecs) results in a best-fit slope that is closer to unity, as observed ($0.96 \pm 0.03$), and increases the dispersion about the relation. The redistribution resulting from this transformation implies that the best-fit parameters and dispersion of the intrinsic $D_A^2 Y_{SZ}-Y_X$ relation, and by extension the $M_{500}^{SZ}-M_{500}^{Y_X}$ relation derived from these quantities, will be sample-dependent.
By itself, our study cannot estimate the absolute value of any bias between the Planck mass and the true mass, or exclude any mass dependence of this bias. Such a mass dependence has been suggested by several lensing studies (von der Linden et al. 2014; Hoekstra et al. 2015) although the latest compilation study by Sereno & Ettori (2017) suggests that the apparent mass-dependent bias is actually due to an underlying redshift dependence. All these works are based on subsamples of Planck clusters, with no simple selection criteria, so that selection effects are potentially complex. To make further progress in this field, we require a complete follow-up of a well-characterized Planck cluster sample with mass estimation using various techniques (X-ray proxies, HE masses, and lensing masses). Such data will be available through the CHEX-MATE program (The CHEX-MATE Collaboration 2020).

F.A.-S. acknowledges support from Chandra grant GO3-14131X. G.W.P., M.A., and J.-B.M. acknowledge funding from the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement No. 340519, and from the French Space Agency, CNES. W.F., C.J., A.V., R.K., and L.D. acknowledge support from the Smithsonian Institution and the Chandra HRC Project through NASA contract NAS8-03060. R.J.v.W. acknowledges support from the VIDI research program with project No. 639.042.729, which is financed by the Netherlands Organization for Scientific Research. L.L. acknowledges support from the contracts ASI-INAF Athena 2019-27-HH.0, Attivitàdi Studio per la Comunità Scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare (Accordo Attuativo ASI-INAF n.2017-14-H.0), and INAF Call per Interventi Aggiuntivia Sostegno della Ricerca di Main-stream di INAF and financial contributions from NASA through contracts 80NSSCK0582 and 80NSSC19K0116. Basic research in radio astronomy at the Naval Research Laboratory is supported by 6.1 Base funding. S.B. acknowledges financial support from the agreement ASI-INAF n.2017-14-H.0 and an INFN INDARK grant.

# Appendix A
## Sample Selection Bias, Scatter, and Covariance Tests

### A.1. Simulations

The present clusters were selected from the ESZ catalog, produced from detections on the Planck 10 month maps, while their PSZ2 $Y_{SZ}$ flux values were extracted from the full-mission maps. To quantify any possible residual bias due to the sample selection, we undertook an extensive series of simulations by generating mock ESZ and PSZ2 catalogs.

Starting from a Tinker mass function (Tinker et al. 2008), simulated clusters were injected into the Planck SZ “cosmological” mask region (excluding the Galactic plane and the Magellanic Clouds) of the 10 month maps used for the original ESZ catalog construction. Each simulated object was modeled with the Arnaud et al. (2010) pressure profile. The $Y_{SZ}$ value corresponding to the cluster redshift and mass was drawn from the $D_A^2 Y_{SZ} - M_{500}$ relation, with a bias between the X-ray-calibrated mass and the true mass, $(1-b) = 0.65$, tuned to recover the cluster number counts for the Planck cosmology.

For the $Y_{X}$ value, we considered $Y_{X}$ physical quantities normalized by the $c^*$ factor given by Equation (14), i.e., expressed in the same units as the intrinsic Compton parameter. We note that any selection effect will not depend on the mean

ratio between the true $Y_{X}$ and $D_A^2 Y_{SZ}$ values. The slope of the $D_A^2 Y_{SZ} - Y_{X}$ relation is expected to be close to unity, and we are interested in possible slope biases induced by selection and measurements. For simplicity, we thus assume that the normalization and slope of the $Y_{X}-M_{500}$ and $D_A^2 Y_{SZ} - M_{500}$ scaling relations are the same.

However, we allow for different scatters and a possible covariance between $Y_{X}$ and $D_A^2 Y_{SZ}$ variations from the relation. We can thus write

$$Y_M = AE(z)^{2/3}M^B \\ P(D_A^2 Y_{SZ}, Y_X|Y_M) = N(Y_M, V_\sigma) \quad (A1)$$

where $Y_M$ is the latent $Y$ value at a given mass (i.e., the values of $Y_X$ and $D_A^2 Y_{SZ}$ obtained if there were no scatter), and $D_A^2 Y_{SZ}$ and $Y_X$ are the true values. We assume a Gaussian lognormal correlated distribution for $Y_{SZ}$ and $Y_X$ at fixed mass, $N(Y_M, V_\sigma)$, with a covariance matrix, $V_\sigma$:

$$\begin{pmatrix} \sigma_{\log, Y_{SZ}}^2 & r \sigma_{\log, Y_{SZ}} \sigma_{\log, Y_X} \\ r \sigma_{\log, Y_{SZ}} \sigma_{\log, Y_X} & \sigma_{\log, Y_X}^2 \end{pmatrix}$$

where $\sigma_{\log, Y_{SZ}}$ and $\sigma_{\log, Y_X}$ are the scatter in the log-log plane of the Compton parameter and $Y_X$ at fixed mass, respectively. The intrinsic $Y_X$ and $D_A^2 Y_{SZ}$ quantities were converted to flux using the angular distance to each cluster. In the absence of bias, the observed $Y_{SZ} - D_A^{-2} Y_X$ relation should be the identity, i.e., it should have a slope and normalization equal to unity. In the above, all quantities refer to values computed within a sphere of radius $R_{500}$.

SZ detections were then obtained by applying the MMF algorithm (Melin et al. 2006) to the 10 month maps. The mock ESZ catalogs were constructed from the application of an $S/N > 6$ threshold to these detections. These mock ESZ clusters were then reinjected into the full-mission maps. The MMF algorithm was then run on these full-mission maps, and the mock PSZ2 value at the true position was extracted for the mock ESZ clusters. We extracted the mock PSZ2 value both at the true size and at the X-ray size. The measurement errors were estimated from the maps as described in Melin et al. (2006). The simulations are therefore fully representative of the observations. We ran 10 ESZ simulations, yielding a total number of 1657 clusters and allowing a very precise estimate of any bias on the slope (the precision is better than ~0.5%). For simplicity we assume that the statistical errors on $Y_X$ are negligible.

In the following, we detail a step-by-step analysis of the simulated relations and the interdependence of Malmquist bias, intrinsic scatter, and covariance between the quantities. We focus on the $Y_{SZ} - D_A^{-2} Y_X$ relation between the SZ flux and its X-ray equivalent. This is the most fundamental observed relation, on which selection effects are expected to be most visible (particularly at the low end). The left-hand panel of Figure 9 shows the effect of Malmquist bias on the observed $Y_{SZ}$ values as compared to the latent $D_A^{-2} Y_X$ value, while the other panels show the corresponding $Y_{SZ} - D_A^{-2} Y_X$ relation and the effect of $Y_X$ intrinsic scatter and covariance with $D_A^2 Y_{SZ}$. Unless otherwise stated the fit results were obtained with the BCES method.

### A.1.1. Malmquist Bias due to Planck Noise and Intrinsic Scatter in the $D_A^2 Y_{SZ} - M_{500}$ Relation

We first assume that there is no intrinsic scatter in the $D_A^2 Y_{SZ} - M_{500}$ relation. The corresponding observed $Y_{SZ}$ value, extracted at
the true size, is shown in the top left panel of Figure 9 as a function of $D_A^{-2}Y_M$. The scatter in the Y-axis direction is entirely due to the observational uncertainties in the SZ measurements. Object selection was performed on the 10 month (ESZ) map and the $Y_{SZ}$ signal was extracted from the full (29 month) mission maps. The noise between the ESZ and PSZ2 maps is slightly correlated, so we estimate here the residual effect of the ESZ selection on the flux estimation in the PSZ2 map. The effect on the slope is very small (slope $0.994 \pm 0.003$, implying a negligible impact.

We now include intrinsic scatter in the $D_A^2 Y_{SZ} - M_{500}$ relation, assuming $\sigma_{\log Y_{SZ}} = 0.05$ (Kay et al. 2012; Le Brun et al. 2017). The addition of intrinsic scatter increases the dispersion in the Y-direction, and has the effect of changing the sample selection. The resulting relation is shown in the bottom left panel of Figure 9. Objects that are newly detected as compared to the previous case (about 10%) are plotted in red, and objects that are lost (no longer detected) are plotted in green (7%). As expected, lost objects are on average less intrinsically bright than newly detected objects, especially at low flux, as more upscattered systems pass the threshold than downscattered systems. The intrinsic scatter then pushes the slope slightly away from unity, to $0.98 \pm 0.005$ (a 2% effect). This is however still far from the observed value of $0.89 \pm 0.01$.

### A.1.2. Using the $Y_X$ Proxy and the Eddington Bias

We now consider $D_A^{-2}Y_X$ as a covariate, the observable we consider in this study. We assumed $\sigma_{\log Y_X} = 0.07$ following numerical simulations (Planelles et al. 2014; Le Brun et al. 2017; Truong et al. 2018). The corresponding $Y_{SZ} - D_A^{-2}Y_X$ relations are shown in the middle panel of Figure 9, for the two cases described in the previous section. $Y_X$ is a scattered estimate of $Y_M$, so the effect is to introduce scatter in the X-axis direction. The extraction at the X-ray size rather than at the true size changes slightly the $Y_{SZ}$ value, but the effect is very small. We also note that the slope is decreased, as compared to that of the $Y_{SZ} - Y_M$ relation, by about $\Delta B = -0.02$, with a slope of $0.972 \pm 0.005$ for $\sigma_{\log Y_{SZ}} = 0$ and of $0.962 \pm 0.006$ for $\sigma_{\log Y_{SZ}} = 0.05$.

This slope decrease is reminiscent of Eddington-like bias, as described extensively by Sereno & Ettori (2015). This bias is due to the intrinsic scatter of the covariate (here $Y_X$) with respect to the latent quantity ($Y_M$) when it has a nonuniform distribution. However this does not fully explain the effect, which we fail to fully understand. Indeed, the same slope is found using the LIRA regression code, designed to correct for this effect: $0.95 \pm 0.02$, for the nominal case with intrinsic scatter in $Y_{SZ}$. The same agreement between LIRA and BCES results was noted for the observations (Table 3).

We next added covariance between the $Y_X$ and $Y_{SZ}$ deviations, assuming a correlation coefficient $r=0.4$ (e.g., Farahi et al. 2019; Nagarajan et al. 2019). The bottom right panel of Figure 9 shows the resulting relation, indicating that addition of covariance pushes the slope back toward unity, to $0.973 \pm 0.005$. This is due to the correlation between the deviations of the quantities in the two axes. Increasing the correlation coefficient from $r=0.4$ to $r=0.6$ increases the push toward unity with the slope changing to $0.982 \pm 0.005$. For a perfect correlation, $Y_X$ and $Y_{SZ}$ would simply move along the line, canceling the Malmquist bias due to intrinsic scatter.

### A.1.3. Extreme Case: Large Value of Intrinsic Scatter in the $M_{500}-Y_{SZ}$ Relation

As a final test, we doubled the intrinsic scatter in the $M_{500}-Y_{SZ}$ relation, so that $\sigma_{\log Y_{SZ}} = 0.10$. If $r=0$, the increased dispersion in the Y-direction again affects the sample selection, but more strongly than in the bottom left panel of Figure 10, driving the slope to $0.930 \pm 0.008$. However, as in the nominal case, inclusion of covariance again redresses the slope toward unity: for $r=0.4$, the slope becomes $0.955 \pm 0.007$, and for $r=0.6$, the slope is $0.967 \pm 0.006$.

The general conclusion is that biases due to selection effects and intrinsic scatter cannot explain the observed slope of the $Y_{SZ}-D_A^{-2}Y_X$ flux relation. As we used PSZ2 SZ values, largely independent of the detection values, the residual Malmquist bias is dominated by the intrinsic scatter between $D_A^2 Y_{SZ}$ and the mass. The intrinsic scatter in $Y_X$ is found to slightly further reduce the slope, but covariance between $Y_X$ and $D_A^2 Y_{SZ}$ redresses the slope toward unity. For typical scatter and covariances derived from simulations, the net effect is a decrease of the slope to $B \sim 0.97$, far from the observed value of $0.89 \pm 0.01$. The combination of factors that gives the most deviant slope from unity requires unrealistically large intrinsic scatter in the $D_A^2 Y_{SZ} - M_{500}$ relation and zero covariance between $Y_X$ and $D_A^2 Y_{SZ}$, and even when doubling the $D_A^2 Y_{SZ}$ scatter the slope does not agree with the data.

### A.2. Further Tests

We further checked that the Malmquist bias due to Planck noise fluctuations is indeed negligible in the data when using PSZ2 values for the ESZ clusters, by examining the correlation between the PSZ2 $Y_{SZ}/Y_X$ and the S/N of the ESZ detection. The $Y_{SZ}$ measured in the ESZ maps is affected by Malmquist bias, as clusters with a signal that is boosted by positive noise have a higher probability of being detected, particularly close to the detection threshold. The signature of this bias is a correlation between the ratio of the observed signal to the true signal, $Y_{SZ}/Y_{SZ, true}$, and the S/N of the ESZ detection. If the PSZ2 measurements were fully independent, this correlation would disappear. We do not have access to $Y_{SZ, true}$, but we can use $Y_X$ as a proxy as it is not affected by Planck noise in the SZ detection. A negligible dependence of $Y_{SZ}/Y_X$ on S/N would imply negligible Malmquist bias due to Planck noise.

However, a difficulty is that more massive clusters are easier to detect, so that the S/N is correlated to mass as shown on the left panel of Figure 10. To disentangle a possible intrinsic mass dependence of $Y_{SZ}/Y_X$ from a dependence on the S/N due to the Malmquist bias, we selected a subsample of clusters of nearly the same mass in a very restricted mass range of $M_{500}^Y = (8.5 \pm 1.5) \times 10^{14} M_{\odot}$. The mean mass of the subsample was chosen to maximize the S/N leverage, while $\Delta M$ is small enough to ensure that there is no residual correlation between mass and S/N. The clusters in this subsample are identified by the red points in the left panel of Figure 10. $Y_{SZ}/Y_X$ for the full sample is significantly correlated with the S/N and increases at low S/N (central panel of Figure 10). However the correlation with S/N disappears for the subsample at a given mass range ($r = -0.015, p = 0.92$). The same is found for the mass ratio ($r = -0.025, p = 0.87$). The correlation coefficients for the full sample and subsample are also given in the plots of Figure 10 and in Table 4.

This test shows that the Malmquist bias due to Planck noise fluctuations is indeed negligible when considering PSZ2 $Y_{SZ}$
values. Note that this does imply a negligible Malmquist bias due to the intrinsic scatter of $Y_{SZ}$ with mass. Clusters detected because their $Y_{SZ}$ are scattered up at a given mass may also have a higher $Y_X$ value due to covariance between $Y_{SZ}$ and $Y_X$ deviations at a given mass, so that the ratio $Y_{SZ}/Y_X$ remains weakly dependent on S/N.

## Appendix B
### Testing Mass Dependence on Chandra-Planck Center Offset

The uncertainty on the Planck cluster center is determined by the spatial resolution of the Planck instruments. Here, we compare the cluster centers determined by the Planck detection algorithm (~arcminutes) to the precise (~arcsecond) Chandra centroids. The mean offset between the Chandra and Planck centers is 1.7, similar to the result obtained by Planck Collaboration IV (2013) of 1.5, and in agreement with expectations from the Planck sky simulations (Planck Collaboration VIII 2011, see their Section 6.1).

We investigate a possible dependence of the SZ-derived to X-ray-derived mass ratio (see Figure 11) on the offset between the Planck and Chandra centers. The relation between the SZ-derived to X-ray-derived mass ratio and the offset between the Planck and Chandra positions ($\Delta R$) is given by

$$M_{500}^{SZ}/M_{500}^{Y_X} = 10^A \times \Delta R^B. \quad (B1)$$

Our results are presented in Table 5. We find no evidence of a correlation between the two quantities (under a null hypothesis of no correlation, we obtained a p-value of 0.46), suggesting that the large Planck position uncertainty (~2') is not driving the differences in cluster mass determinations.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Ratio between SZ-derived and X-ray-derived masses plotted against the projected distances between the cluster centers as measured using Planck and Chandra data. The black solid line corresponds to the power-law best fit, while the red solid line corresponds to equal SZ-derived and X-ray-derived masses. The dashed line corresponds to the 1σ confidence range. The dotted line corresponds to the 1σ envelope on the best-fitting line.</figcaption></figure>

<table><caption>**Table 5**
Best-fit Parameters for the Ratio between SZ and X-Ray Masses</caption><thead><tr><th>Relation</th><th>A</th><th>σ<sub>A</sub></th><th>B</th><th>σ<sub>B</sub></th></tr></thead><tbody><tr><td>Ratio versus ΔR</td><td>-0.038</td><td>0.006</td><td>-0.027</td><td>0.025</td></tr></tbody></table>

**Note.** Columns list the best-fit parameters and their uncertainties for the power law given by Equation (B1).

## Appendix C
### Difference in Slope of the $Y_{SZ}-Y_X$ Relation Expressed in Apparent Flux (Units of Square Arcminutes) and Intrinsic Properties (Units of Square Megaparsecs)

As noted throughout the text, we observe a statistically significant difference in the slope of the $Y_{SZ}-Y_X$ relation when it is expressed in apparent flux (units of square arcminutes) compared to intrinsic properties (units of square megaparsecs), as clearly seen also in Figure 12. When clusters are color-coded according to their gas temperatures, we clearly see on the right panel of Figure 12 that the clusters are ordered by color, as expected when the $Y_{SZ}-Y_X$ relation is expressed in intrinsic properties (units of square megaparsecs). On the other hand, on the left panel of Figure 12 the colors are more evenly spread throughout the $Y_{SZ}-Y_X$ relation. A reordering of the data points occurs when the $Y_{SZ}-Y_X$ relation is expressed in intrinsic properties instead of apparent flux because this change depends on the underlying cosmology and cluster redshift, the result being a different multiplying factor for each cluster. In this appendix we show that this change from apparent flux to intrinsic properties results in a best-fit $Y_{SZ}-Y_X$ relation slope that is closer to unity when expressed in intrinsic properties (units of square megaparsecs).

### C.1. Analytical Demonstration of the Change in Slope and Dispersion

For simplicity, the linear regression presented here is undertaken using the least squares method. For *n* pairs of data (x, y), the method of least squares may be used to write a linear relationship between *x* and *y*. The least squares regression line is the locus that minimizes the sum of the squares of the vertical deviation from each data point from the best-fit linear relation.

The least squares regression line for a set of *n* data points is given by the equation of a line:

$$y = a + bx, \quad (C1)$$

where *a*, *b*, and the variance, $\sigma^2$, about the best fit are given by ($\sum_{i=1}^n \rightarrow \Sigma$)

$$b = \frac{n \sum x_i y_i - (\sum x_i)(\sum y_i)}{n \sum x_i^2 - (\sum x_i)^2} \quad (C2)$$

$$a = \frac{\sum y_i - b \sum x_i}{n} \quad (C3)$$

$$\sigma^2 = \frac{\sum (y_i - a - b x_i)^2}{n - 1}. \quad (C4)$$

Assuming that the data are now best described by a power-law relation given by $y = ax^b$, one can linearize this relation by applying a logarithmic function to both sides of this equality:

$$\log(y) = \log(a) + b \log(x) \rightarrow z = a_0 + a_1 w. \quad (C5)$$
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Observed change in slope in the $Y_{SZ}-Y_X$ relation caused by the change from apparent flux to intrinsic properties. The points are color-coded according to the values of the gas temperature. As we can see, while the colors are separated in the $Y_{SZ}-Y_X$ relation in intrinsic properties (units of square megaparsecs), they are mixed when the same relation is expressed in apparent flux (units of square arcminutes), illustrating the redistribution that happens when the $Y_{SZ}-Y_X$ relation is expressed in different physical quantities. This redistribution ultimately changes the slope of the $Y_{SZ}-Y_X$ relation. Left: The best-fit relation is given by $Y_{SZ} \propto (Y_X)^{0.89}$ in apparent flux (units of square arcminutes). Right: Same as left panel, but for $Y_{SZ} \propto (Y_X)^{0.96}$ in intrinsic properties (units of square megaparsecs).</figcaption></figure>

Now, let us perform the following linear transformation on the data (x, y):

$$y' = \zeta y: y'_{i} = \zeta_{i} y_{i} \text{ and } x' = \zeta x: x'_{i} = \zeta_{i} x_{i}, \quad (C6)$$

where the coefficient $\zeta_i$ is unique for each pair of data $(x_i, y_i)$. If each data point $(x_i, y_i)$ has a dispersion about the best-fit relation given by $1 + \epsilon_i$, such that $y_i = ax_i^b(1 + \epsilon_i)$, the linear transformation leads to $y_i' = a\zeta_i x_i^b(1 + \epsilon_i) = a\zeta_i^{1-b}x_i^b(1 + \epsilon_i)$, which can be linearized as

$$\begin{aligned} \log(y_i') &= \log(a) + (1-b)\log(\zeta_i) \\ &\quad + \log(1+\epsilon_i) + b \log(x_i') \rightarrow z_i' \\ &= a_0 + (1-b)a_{1,i} + a_{2,i} + bw_i'. \end{aligned} \quad (C7)$$

Using the linearized equation presented in (C7), Equation (C3) takes the following form:

where the term multiplying $(1-b)$ is the slope of the best fit of $\zeta(x') \equiv \zeta(\zeta x)$ in logarithmic space (compare Equation (C10) with Equation (C2)), which will be named $\beta$. The last term is the slope of the best fit of the $(1+\epsilon_i)$-$\zeta x$ relation in logarithmic space. Assuming there is no correlation between the dispersion about the original best-fit relation and the transformed x-coordinate ($\zeta x$), the last term of Equation (C9) vanishes, so that it finally becomes

$$b' = b + (1-b) \frac{n \sum \log(\zeta_i) \log(\zeta_i x_i) - \sum \log(\zeta_i) \sum \log(\zeta_i x_i)}{n \sum \log(\zeta_i x_i)^2 - (\sum \log(\zeta_i x_i))^2}. \quad (C10)$$

Following the same arguments, we can also express (for completeness) the new linear coefficient $a'$ as

$$a' = \log(a) + (1-b) \frac{\sum \log(\zeta_i) - \beta \sum \log(\zeta_i x_i)}{n}, \quad (C11)$$

$$\begin{aligned} b' &= \frac{n \sum w_i' z_i' - (\sum w_i')(\sum z_i')}{n \sum w_i'^2 - (\sum w_i')^2} \\ &= \frac{n \sum (a_0 w_i' + (1-b)a_{1,i}w_i' + a_{2,i}w_i' + bw_i'^2) - (\sum w_i')(\sum (a_0 + (1-b)a_{1,i} + a_{2,i} + bw_i'))}{n \sum w_i'^2 - (\sum w_i')^2} \\ &= \frac{n \sum a_0 w_i' + n \sum (1-b)a_{1,i}w_i' + n \sum a_{2,i}w_i' + n \sum bw_i'^2 - n \sum a_0 w_i' - \sum (1-b)a_{1,i} \sum w_i' - \sum a_{2,i} \sum w_i' - b(\sum w_i')^2}{n \sum w_i'^2 - (\sum w_i')^2} \\ &= b + (1-b) \frac{n \sum a_{1,i}w_i' - \sum a_{1,i} \sum w_i'}{n \sum w_i'^2 - (\sum w_i')^2} + \frac{n \sum a_{2,i}w_i' - \sum a_{2,i} \sum w_i'}{n \sum w_i'^2 - (\sum w_i')^2}, \end{aligned} \quad (C8)$$

which becomes

$$b' = b + (1-b) \frac{n \sum \log(\zeta_i) \log(\zeta_i x_i) - \sum \log(\zeta_i) \sum \log(\zeta_i x_i)}{n \sum \log(\zeta_i x_i)^2 - (\sum \log(\zeta_i x_i))^2} + \frac{n \sum \log(1+\epsilon_i) \log(\zeta_i x_i) - \sum \log(1+\epsilon_i) \sum \log(\zeta_i x_i)}{n \sum \log(\zeta_i x_i)^2 - (\sum \log(\zeta_i x_i))^2}, \quad (C9)$$

where the term multiplying $(1-b)$ is the linear coefficient of the best fit of $\zeta(x') \equiv \zeta(\zeta x)$ in logarithmic space (compare Equation (C11) with Equation (C3)), which will be named $\alpha$.

We can finally express the new dispersion $\sigma'$, following again steps similar to those leading to Equation (C10), by
<figure><img src="image_16.png" /><figcaption>Figure 13. Visualization of the $\zeta-\zeta x$ relation in logarithmic space. Top left: Illustration of boundaries for the pairs $(\zeta_i, \zeta_i x_i)$ when $\zeta_i \in [\zeta_{\text{inf}}, \zeta_{\text{sup}}]$ and $x_i \in [x_{\text{inf}}, x_{\text{sup}}]$. Top right: Illustration of the slope of the best power-law fit when $x_{\text{inf}} \sim x_{\text{sup}}$ and $\zeta_{\text{inf}} \ll \zeta_{\text{sup}}$. Bottom left: Same as top right panel, but for $x_{\text{inf}} \ll x_{\text{sup}}$ and $\zeta_{\text{inf}} \sim \zeta_{\text{sup}}$. Bottom right: Same as top right panel, but for $x_{\text{inf}}/x_{\text{sup}} = \zeta_{\text{inf}}/\zeta_{\text{sup}}$.</figcaption></figure>

$$ \sigma'^2 = \frac{\sum \log(1 + \epsilon_i)^2}{n-1} + (1-b)^2 \frac{\sum \log(\zeta_i - \alpha - \beta \log(\zeta_i x_i))^2}{n-1}, \quad (\text{C12}) $$

where the first term is the variance of the starting relation, $\sigma^2$, and the term multiplying $(1-b)^2$ is the variance of $\zeta(x') \equiv \zeta(\zeta x)$ with respect to its best fit in logarithmic space (compare Equation (C12) with Equation (C4)), which will be named $\gamma^2$.

We can now express Equations (C10), (C11), and (C12) simply as

$$ \begin{aligned} a' &= \log(a) + \alpha (1 - b), \\ b' &= b + \beta (1 - b), \\ \sigma'^2 &= \sigma^2 + \gamma^2 (1 - b)^2. \end{aligned} \quad (\text{C13}) $$

From Figure 13, we can easily see that $0 \le \beta \le 1$. This inequality, in conjunction with the set of Equation (C13), leads to the immediate conclusion that for $b \ne 1$ the linear transformation we have just presented will result in a best-fit relation that is closer to unity, while increasing the dispersion of the data around this best fit.

Now let us assume that both $\zeta$ and $x$ are random variables uniformly distributed in the intervals $[\zeta_{\text{inf}}, \zeta_{\text{sup}}]$ and $[x_{\text{inf}}, x_{\text{sup}}]$ (see top left panel of Figure 13). If $x_{\text{inf}} \sim x_{\text{sup}}$ and $\zeta_{\text{inf}} \ll \zeta_{\text{sup}}$, the slope of the best fit of $\zeta(\zeta x)$ is $\sim 1$ (see top right panel of

Figure 13), and Equation (C10) becomes $b' \approx 1$. If $x_{\text{inf}} \ll x_{\text{sup}}$ and $\zeta_{\text{inf}} \sim \zeta_{\text{sup}}$, the slope of the best fit of $\zeta(\zeta x)$ is $\sim 0$ (see bottom left panel of Figure 13), and Equation (C10) becomes $b' \approx b$. If $x_{\text{inf}}/x_{\text{sup}} = \zeta_{\text{inf}}/\zeta_{\text{sup}}$, the slope of the best fit of $\zeta(\zeta x)$ is $1/2$ (see bottom right panel of Figure 13), and Equation (C10) becomes $b' = b + (1-b)/2$, which means that after random distributions of $x$ and $y$ (with the same factor, $\zeta_i$, multiplying both coordinates for each data point) the best-fit slope is closer to unity, decreasing the difference by half. To test this analytical prediction we created mock representations of the $Y_{SZ}-Y_X$ relation, which are presented in the next section.

## C.2. Mock Representations of the $Y_{SZ}-Y_X$ Relation

To test the analytical prediction presented in the previous section, we created mock representations of the $Y_{SZ}-Y_X$ relation, using 200 data points. Initially, we populated the $Y_{SZ}-Y_X$ relation with values in a range similar to that presented in Figure 12 to mimic the $Y_{SZ}-Y_X$ relation expressed in apparent flux (units of square arcminutes). We then randomly redistributed the data points using a random variable uniformly distributed in the [0, 0.1] range to mimic the values presented in the $Y_{SZ}-Y_X$ relation expressed in intrinsic properties (units of square megaparsecs). The mock representations of the $Y_{SZ}-Y_X$ relation are presented in Figure 14. As we can see, the best-fit results agree perfectly with the analytical prediction ($b' = b + (1-b)/2$) presented in Section C.1 (see Table 6).
<figure><img src="image_17.png" /><figcaption>Figure 14. Toy model of the change in slope caused by the change of units from square arcminutes to square megaparsecs. The points are color-coded according to the mock values of $Y_X$ in units of square megaparsecs, which should be closely related to a temperature dependence. As we can see, while the colors are separated in the $Y_{SZ}-Y_X$ relation in units of square megaparsecs, they are mixed when the same relation is plotted in units of square arcminutes, illustrating the redistribution that happens when the values are presented in other units. This redistribution ultimately changes the slope of the best fit in this particular set of examples as given by $b' = b + (1 - b)/2$. Top left: Creating a distribution that follows $Y_{SZ} \propto (Y_X)^1$ in units of square arcminutes, one gets the same slope in units of square megaparsecs, as expected from $b' = b + (1 - b)/2$. Top right: Same as top left panel, but for $Y_{SZ} \propto (Y_X)^{0.90}$ in units of square arcminutes leading to $Y_{SZ} \propto (Y_X)^{0.95}$ in units of square megaparsecs. Bottom left: Same as top left panel, but for $Y_{SZ} \propto (Y_X)^{0.80}$ in units of square arcminutes leading to $Y_{SZ} \propto (Y_X)^{0.90}$ in units of square megaparsecs. Bottom right: Same as top left panel, but for $Y_{SZ} \propto (Y_X)^{0.70}$ in units of square arcminutes leading to $Y_{SZ} \propto (Y_X)^{0.86}$ in units of square megaparsecs.</figcaption></figure>

### C.3. Predicted versus Observed Change in Slope for Our Sample

Finally, we investigate the predicted change in slope for our sample, using the factors ($\zeta_i$) that were used for the conversion from apparent flux (units of square arcminutes) to intrinsic properties (units of square megaparsecs). Using the notation of Section C.1 and starting from the $Y_{SZ}-Y_X$ relation expressed in apparent flux (units of square arcminutes), we have

$$1. x_i = Y_{X,arcmin^2}^i \text{ and } y_i = Y_{SZ,arcmin^2}^i;$$

$$2. \zeta_i = Y_{SZ,Mpc^2}^i / Y_{SZ,arcmin^2}^i = Y_{X,Mpc^2}^i / Y_{X,arcmin^2}^i \text{ and } \zeta_i x_i = Y_{X,Mpc^2}^i; \text{ and}$$

$$3. b = \text{slope of the best fit of the } Y_{SZ}-Y_X \text{ relation in apparent flux (units of square arcminutes)} \text{ and } b' = \text{slope of the best fit of the } Y_{SZ}-Y_X \text{ relation in intrinsic properties (units of square megaparsecs).}$$

Using the formalism presented in Section C.1 and the variables defined above, the predicted slope of the best fit of the $Y_{SZ}-Y_X$ relation when it is changed from apparent flux (units of square arcminutes) to intrinsic properties (units of square megaparsecs) is
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Top: Observed change in slope in the $Y_{SZ}-Y_X$ relation caused by the change of units from square arcminutes to square megaparsecs. The points are color-coded according to the cluster redshift. As we can see, while the colors are separated in the $Y_{SZ}-Y_X$ relation in units of square arcminutes, they are mixed when the same relation is plotted in units of square megaparsecs, showing a contrary effect compared to when the colors are coded according to the cluster gas temperature (when the points are color-coded according to the gas temperature, the colors are separated in units of square megaparsecs—Figure 12). This illustrates the redistribution that happens when the values are presented in different units. This redistribution ultimately changes the slope of the $Y_{SZ}-Y_X$ relation. Top left: The best-fit relation is given by $Y_{SZ} \propto (Y_X)^{0.89}$ in units of square arcminutes. Top right: Same as top left panel, but for $Y_{SZ} \propto (Y_X)^{0.96}$ in units of square megaparsecs. Bottom: Observed $\zeta-\zeta_X$ relation in logarithmic space for our ESZ sample. As in the top panels, the points are color-coded according to the cluster redshift. Bottom left: Observed $\zeta-\zeta_X$ relation used for predicting the change of slope of the $Y_{SZ}-Y_X$ relation when units are changed from square megaparsecs to square arcminutes. The best-fit slope is given by $\beta = 0.90 \pm 0.07$. Bottom right: Same as bottom left panel, but for changing of units from square arcminutes to square megaparsecs. The best-fit slope is given by $\beta = 0.85 \pm 0.10$.</figcaption></figure>

given by Equation (C10), where the term multiplying $(1-b)$ is the slope of the best fit of the $(Y_{SZ}, \text{Mpc}^2 / Y_{SZ, \text{arcmin}^2}) - Y_X, \text{Mpc}^2$ relation ($\zeta-\zeta_X$) in logarithmic space, which turns out to be 0.85 (see bottom panel of Figure 15). This leads to a predicted slope of the best fit of the $Y_{SZ}-Y_X$ relation in intrinsic properties of $b' = 0.89 + (1-0.89) \times 0.85 = 0.98$, which is in agreement with the observed value of $0.96 \pm 0.03$. Now, if we perform the same calculation to compute the predicted value of the slope of $Y_{SZ}-Y_X$ in apparent flux starting from this relation in intrinsic properties (the slope of the best fit of $\zeta-\zeta_X$ is 0.90—see bottom

panel of Figure 15), we obtain $b' = 0.96 + (1-0.96) \times 0.90 = 1.00$, which is in total disagreement with the observed value of $0.89 \pm 0.01$.

The question that then arises naturally is why going from apparent flux to intrinsic properties results in a slope that is closer to unity, and yet the reverse is not the case.

Based on the analytical formalism presented in Section C.1, we know that any reordering of the data will result in a slope that is closer to unity. However, the analytical formalism does not indicate any mathematical property that allows us to know
<table><caption>**Table 6**
Best-fit Slope for $Y_{SZ}-Y_X$ Expressed in Apparent Flux (Units of Square Arcminutes) and Intrinsic Properties (Square Megaparsecs) for Our Mock Relations</caption><thead><tr><th>B<sub>arcmin</sub><sup>2</sup></th><th>σ<sub>B</sub><sub>arcmin</sub><sup>2</sup></th><th>B<sub>Mpc</sub><sup>2</sup></th><th>σ<sub>B</sub><sub>Mpc</sub><sup>2</sup></th><th>Predicted</th></tr></thead><tbody><tr><td>1.00</td><td>0.01</td><td>1.00</td><td>0.01</td><td>1.00</td></tr><tr><td>0.90</td><td>0.01</td><td>0.95</td><td>0.01</td><td>0.95</td></tr><tr><td>0.80</td><td>0.01</td><td>0.90</td><td>0.01</td><td>0.90</td></tr><tr><td>0.70</td><td>0.01</td><td>0.86</td><td>0.02</td><td>0.85</td></tr></tbody></table>

**Note.** Columns list best-fit slopes and their uncertainties for the power-law fit to the $Y_{SZ}-Y_X$ relation expressed in apparent flux (units of square arcminutes) and intrinsic properties (units of square megaparsecs) and the predicted value given by $B_{Mpc}^2 = B_{arcmin}^2 + (1 - B_{arcmin}^2)/2$.

the direction in which this transformation results in a slope that is closer to unity.

The answer to this question seems to be given by understanding what is first observed—that is, by understanding what is the starting point. The observed integrated Compton parameter ($Y_{SZ}$) is measured from the apparent flux (square arcminutes) from the Planck data. Similarly, the X-ray equivalent of the integrated Compton parameter, $Y_X$, which is given by the product of the gas mass and gas temperature, has a cosmological dependence on the gas mass, which is originally measured from the observed emission measure profile, and which is therefore also dependent on the apparent flux. It is only after calculations involving the underlying cosmology and cluster redshift that $Y_X$ is presented in intrinsic properties ($M_\odot$ keV or $Mpc^2$). Based on this reasoning, we see that the starting point is given by measurements that involve the apparent flux. Therefore, going from apparent flux (units of square arcminutes) to intrinsic properties (units of square megaparsecs) should result in a best-fit slope that is closer to unity according to our analytical formalism presented in Section C.1. This is indeed what we observe.

## ORCID iDs

Felipe Andrade-Santos https://orcid.org/0000-0002-8144-9285

William R. Forman https://orcid.org/0000-0002-9478-1682

Etienne Pointecouteau https://orcid.org/0000-0001-8488-3645

Iacopo Bartalucci https://orcid.org/0000-0001-7703-9040

Alexey Vikhlinin https://orcid.org/0000-0001-8121-0234

Pasquale Mazzotta https://orcid.org/0000-0002-5411-1748

Stefano Borgani https://orcid.org/0000-0001-6151-6439

Lorenzo Lovisari https://orcid.org/0000-0002-3754-2415

Reinout J. van Weeren https://orcid.org/0000-0002-0587-1660

Ralph P. Kraft https://orcid.org/0000-0002-0765-0511

Simona Giacintucci https://orcid.org/0000-0002-1634-9886

## References

Adam, R., Arnaud, M., Bartalucci, I., et al. 2017, A&A, 606, A64

Akritas, M. G., & Bershady, M. A. 1996, ApJ, 470, 706

Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409

Anders, E., & Grevesse, N. 1989, GeCoA, 53, 197

Andrade-Santos, F., Jones, C., Forman, W. R., et al. 2017, ApJ, 843, 76

Andrade-Santos, F., Lima Neto, G. B., & Laganá, T. F. 2012, ApJ, 746, 139

Andrade-Santos, F., Nulsen, P. E. J., Kraft, R. P., et al. 2013, ApJ, 766, 107

Angulo, R. E., Springel, V., White, S. D. M., et al. 2012, MNRAS, 426, 2046

Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92

Bartalucci, I., Arnaud, M., Pratt, G. W., et al. 2017, A&A, 608, A88

Benson, B. A., de Haan, T., Dudley, J. P., et al. 2013, ApJ, 763, 147

Böhringer, H., Pratt, G. W., Arnaud, M., et al. 2010, A&A, 514, A32

Buote, D. A., & Tsai, J. C. 1996, ApJ, 458, 27

Cavaliere, A., & Fusco-Femiano, R. 1976, A&A, 49, 137

da Silva, A. C., Kay, S. T., Liddle, A. R., & Thomas, P. A. 2004, MNRAS, 348, 1401

Dietrich, J. P., Bocquet, S., Schrabback, T., et al. 2019, MNRAS, 483, 2871

Farahi, A., Mulroy, S. L., Evrard, A. E., et al. 2019, NatCo, 10, 2504

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Forman, W., & Jones, C. 1982, ARA&A, 20, 547

Giles, P. A., Maughan, B. J., Pacaud, F., et al. 2016, A&A, 592, A3

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Herranz, D., Sanz, J. L., Hobson, M. P., et al. 2002, MNRAS, 336, 1057

Hoekstra, H., Herbonnet, R., Muzzin, A., et al. 2015, MNRAS, 449, 685

Itoh, N., Kohyama, Y., & Nozawa, S. 1998, ApJ, 502, 7

Jeltema, T. E., Canizares, C. R., Bautz, M. W., & Buote, D. A. 2005, ApJ, 624, 606

Jones, C., & Forman, W. 1984, ApJ, 276, 38

Jones, C., & Forman, W. 1999, ApJ, 511, 65

Jones, C., Mandel, E., Schwarz, J., et al. 1979, ApJL, 234, L21

Kaiser, N. 1986, MNRAS, 222, 323

Kay, S. T., Peel, M. W., Short, C. J., et al. 2012, MNRAS, 422, 1999

Kelly, B. C. 2007, ApJ, 665, 1489

Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353

Laganá, T. F., Andrade-Santos, F., & Lima Neto, G. B. 2010, A&A, 511, A15

Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2017, MNRAS, 466, 4442

Lovisari, L., Ettori, S., Sereno, M., et al. 2020, A&A, 644, A78

Mantz, A., Allen, S. W., Rapetti, D., & Ebeling, H. 2010, MNRAS, 406, 1759

Melin, J., Bartlett, J. G., & Delabrouille, J. 2006, A&A, 459, 341

Mohr, J. J., Evrard, A. E., Fabricant, D. G., & Geller, M. J. 1995, ApJ, 447, 8

Nagarajan, A., Pacaud, F., Sommer, M., et al. 2019, MNRAS, 488, 1728

Piffaretti, R., Jetzer, P., & Schindler, S. 2003, A&A, 398, 41

Pike, S. R., Kay, S. T., Newton, R. D. A., Thomas, P. A., & Jenkins, A. 2014, MNRAS, 445, 1774

Planck Collaboration Int. V 2013, A&A, 550, A131

Planck Collaboration IV 2013, A&A, 550, A130

Planck Collaboration VIII 2011, A&A, 536, A8

Planck Collaboration IX 2011, A&A, 536, A9

Planck Collaboration X 2011, A&A, 536, A10

Planck Collaboration XI 2011, A&A, 536, A11

Planck Collaboration XX 2014, A&A, 571, A20

Planck Collaboration XXIV 2016, A&A, 594, A24

Planck Collaboration XXIX 2014, A&A, 571, A29

Planck Collaboration XXVII 2016, A&A, 594, A27

Planelles, S., Borgani, S., Fabjan, D., et al. 2014, MNRAS, 438, 195

Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, SSRv, 215, 25

Pratt, G. W., Arnaud, M., Piffaretti, R., et al. 2010, A&A, 511, A85

Pratt, G. W., Croston, J. H., Arnaud, M., & Böhringer, H. 2009, A&A, 498, 361

Rasia, E., Ettori, S., Moscardini, L., et al. 2006, MNRAS, 369, 2013

Reiprich, T. H., & Böhringer, H. 2002, ApJ, 567, 716

Robotham, A. S. G., & Obreschkow, D. 2015, PASA, 32, e033

Rozo, E., Vikhlinin, A., & More, S. 2012, ApJ, 760, 67

Schellenberger, G., & Reiprich, T. H. 2017, MNRAS, 469, 3738

Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30

Sereno, M. 2016, MNRAS, 455, 2149

Sereno, M., & Ettori, S. 2015, MNRAS, 450, 3633

Sereno, M., & Ettori, S. 2017, MNRAS, 468, 3322

Stanek, R., Evrard, A. E., Böhringer, H., Schuecker, P., & Nord, B. 2006, ApJ, 648, 956

Sunyaev, R. A., & Zeldovich, Y. B. 1972, CoASP, 4, 173

The CHEX-MATE Collaboration 2020, arXiv:2010.11972

Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Truong, N., Rasia, E., Mazzotta, P., et al. 2018, MNRAS, 474, 4089

Vikhlinin, A., Burenin, R. A., Ebeling, H., et al. 2009b, ApJ, 692, 1033

Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691

Vikhlinin, A., Kravtsov, A. V., Burenin, R. A., et al. 2009b, ApJ, 692, 1060

Vikhlinin, A., Markevitch, M., Murray, S. S., et al. 2005, ApJ, 628, 655

Voevodkin, A., & Vikhlinin, A. 2004, ApJ, 601, 610

Voit, G. M., Kay, S. T., & Bryan, G. L. 2005, MNRAS, 364, 909

von der Linden, A., Mantz, A., Allen, S. W., et al. 2014, MNRAS, 443, 1973

Willingale, R., Starling, R. L. C., Beardmore, A. P., Tanvir, N. R., & O'Brien, P. T. 2013, MNRAS, 431, 394