# Measurement of the Relativistic Sunyaev-Zeldovich Correction in RX
J1347.5-1145

Victoria L. Butler Equal contribution<EMAIL_ADDRESS>Rochester Institute of
Technology, 1 Lomb Memorial Drive, Rochester, NY 14623, USA  Richard M. Feder
Equal contribution.  California Institute of Technology Division of Physics,
Math, and Astronomy, 1200 East California Boulevard, Pasadena, CA 91125, USA
Tansu Daylan Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139  Department of
Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544
Adam B. Mantz Kavli Institute for Particle Astrophysics and Cosmology,
Stanford University, 452 Lomita Mall, Stanford, CA 94305, USA  Dale Mercado
Rochester Institute of Technology, 1 Lomb Memorial Drive, Rochester, NY 14623,
USA  Alfredo Montaña 1Consejo Nacional de Ciencia y Tecnología, Av.
Insurgentes Sur 1582, Col. Crédito Constructor, Alcaldía Benito Juárez, C.P.
03940, Ciudad de México, México
2Instituto Nacional de Astrofísica Óptica y Electrónica, Luis Enrique Erro 1,
CP 72840, Tonantzintla, Puebla, México
Stephen K. N. Portillo DIRAC Institute, Department of Astronomy, University of
Washington, 3910 15th Ave. NE, Seattle, WA 98195 USA  Jack Sayers California
Institute of Technology Division of Physics, Math, and Astronomy, 1200 East
California Boulevard, Pasadena, CA 91125, USA  Benjamin J. Vaughan Rochester
Institute of Technology, 1 Lomb Memorial Drive, Rochester, NY 14623, USA
Michael Zemcov Rochester Institute of Technology, 1 Lomb Memorial Drive,
Rochester, NY 14623, USA  Jet Propulsion Laboratory, 4800 Oak Grove Drive,
Pasadena, CA 91109, USA  Adi Zitrin Physics Department, Ben-Gurion University
of the Negev, P.O. Box 653, Be’er-Sheva 84105, Israel

###### Abstract

We present a measurement of the relativistic corrections to the thermal
Sunyaev-Zel’dovich (SZ) effect spectrum, the rSZ effect, toward the massive
galaxy cluster RX J1347.5-1145 by combining sub-mm images from Herschel-SPIRE
with mm-wave Bolocam maps. Our analysis simultaneously models the SZ effect
signal, the population of cosmic infrared background (CIB) galaxies, and
galactic cirrus dust emission in a manner that fully accounts for their
spatial and frequency-dependent correlations. Gravitational lensing of
background galaxies by RX J1347.5-1145 is included in our methodology based on
a mass model derived from HST observations. Utilizing a set of realistic mock
observations, we employ a forward modelling approach that accounts for the
non-Gaussian covariances between observed astrophysical components to
determine the posterior distribution of SZ effect brightness values consistent
with the observed data. We determine a maximum a posteriori (MAP) value of the
average Comptonization parameter of the intra-cluster medium (ICM) within
R2500 to be $\langle y\rangle_{2500}$ $=1.56\times 10^{-4}$, with
corresponding 68 per cent credible interval $[1.42,1.63]\times 10^{-4}$, and a
MAP ICM electron temperature of
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ $=22.4$ keV with 68 per cent
credible interval spanning $[10.4,33.0]$ keV. This is in good agreement with
the pressure-weighted temperature obtained from Chandra X-ray observations,
$\langle\textrm{T}_{\textrm{x,pw}}\rangle_{2500}$ $=17.4\pm 2.3$ keV. We aim
to apply this methodology to comparable existing data for a sample of 39
galaxy clusters, with an estimated uncertainty on the ensemble mean
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ at the $\simeq 1$ keV level,
sufficiently precise to probe ICM physics and to inform X-ray temperature
calibration.

Galaxy clusters — Sunyaev-Zeldovich effect — Intracluster medium — RX
J1347.5-1145 — cluster mass function

††facilities: Herschel, HST, CXO, Planck, CSO††software: astropy (Astropy
Collaboration et al., 2013, 2018), corner (Foreman-Mackey, 2016), SZpack
(Chluba et al., 2012, 2013), CIAO (https://cxc.cfa.harvard.edu/ciao/,
Fruscione et al. 2006, SPIE Proc. 6270, 62701V, D.R. Silvia & R.E. Doxsey,
eds.) PCAT-DE (Feder et _al_. 2021, in prep.)

##  1 Introduction

Galaxy clusters are the largest bound objects in the Universe, and gravity is
the dominant process driving their evolution and setting their overall
physical characteristics. For our analysis, these objects are largely self-
similar in their properties (Kaiser, 1986). However, deviations from self-
similarity can occur from processes outside of simple matter aggregation, such
as major mergers between similar mass halos, feedback from active galactic
nuclei (AGN), and gas sloshing within the cluster potential (Zhao, 2003;
McNamara & Nulsen, 2007; Markevitch & Vikhlinin, 2007). In some cases, such as
cosmological studies relying on precise halo mass estimates, it is important
to quantify and understand the effects of these more complicated processes
(Pratt et al., 2019). Such deviations from self-similarity are particularly
evident in the spatially resolved thermodynamics of the gaseous intra-cluster
medium (ICM), which have typically been studied via the thermal Bremsstrahlung
emission in X-rays (Sarazin, 1986). However, these studies are generally more
difficult and/or impractical at higher redshifts and larger cluster-centric
radii, due to $(1+z)^{4}$ cosmological dimming and the density squared
dependence of the X-ray surface brightness (McDonald et al., 2017; Reiprich et
al., 2013). In addition, the X-ray instruments on board Chandra and XMM-Newton
have limited sensitivity to ICM gas above $\sim 10$ keV, and such temperatures
are expected to be relatively common among the population of massive galaxy
clusters, for example, in the shock-heated gas that results from major mergers
(Markevitch et al., 2002).

A complementary method to study ICM thermodynamics is through the Sunyaev-
Zeldovich (SZ) effect signal, which is due to inverse Compton scattering of
Cosmic Microwave Background (CMB) photons with energetic electrons in the ICM
(Sunyaev & Zeldovich, 1972). The electron population has both random thermal
motions and coherent velocities, the latter due to internal turbulence or the
overall peculiar velocity of the cluster. These two properties of the ICM give
rise to a thermal and a kinematic SZ effect. In addition, the ICM is mildly
relativistic, with an average normalized temperature equal to approximately 1%
of the normalized electron mass. This gives rise to generally mild
temperature-dependent relativistic corrections to the thermal SZ effect
spectrum known as the rSZ effect signal (Wright, 1979; Itoh et al., 1998;
Chluba et al., 2012, 2013). The typical signatures of these signals are shown
in Figure 1. These relativistic corrections decrease the magnitude of the SZ
effect signal at $\lesssim 500$ GHz, while boosting the signal at higher
frequencies.

Figure 1: The thermal SZ effect spectrum with and without relativistic
corrections, along with kSZ effect. The thermal SZ effect spectra were
calculated assuming a total Comptonization parameter $y=10^{-4}$ and four
different temperature values. As temperature increases, the overall impact of
the rSZ effect is to shift the thermal SZ effect spectrum to higher
frequencies. The kSZ effect was determined using the same value of $y$
combined with a peculiar velocity of $+3000$ km s-1 and an electron
temperature of 17 keV.

While X-ray measurements are a more mature probe, with current satellites
providing spectroscopic imaging with orders of magnitude better statistics and
angular resolution than typical SZ effect observations (Fabian et al., 2011),
the SZ effect can play a critical role in several important areas of study.
First, because the SZ effect is a fractional distortion of the CMB, the
surface brightness is independent of redshift, and so it can more uniformly
probe objects across cosmic time (Bleem et al., 2015). In addition, since the
SZ effect signal is due to photon-electron scattering rather than the ion-
electron acceleration that sources the X-ray Bremmstrahlung, it is relatively
brighter in the lower density cluster outskirt regions (Planck Collaboration
et al., 2013). Finally, the rSZ effect signal is sensitive to arbitrarily
large temperatures above the band limitations of facilities such as Chandra
and XMM-Newton.

While there are numerous examples of SZ effect studies pushing to new regimes
at high redshift and/or large cluster-centric radii, relatively little effort
has been directed toward measuring ICM temperatures with the rSZ effect. This
is due to the combination of sensitivity needed to measure the rSZ effect
signal and the myriad spatial and spectral data required to separate
contamination from other astronomical sources. For example, the contrast
between the rSZ effect signal and the canonical tSZ effect signal is maximized
at higher frequencies in the sub-millimeter, where the total extragalactic
signal is dominated by the Cosmic Infrared Background (CIB, Hauser et al.,
1998).

The difficulty in measuring the rSZ effect is evidenced by the relatively
modest constraints obtained to date. The first measurement to indicate a
deviation from the non-relativistic tSZ effect spectrum (at $\simeq 2\sigma$)
was reported by Zemcov et al. (2010) based on SPIRE observations of the Bullet
cluster (see also Prokhorov et al., 2011; Prokhorov & Colafrancesco, 2012).
Subsequent observations using the Z-Spec spectrometer obtained a similar
detection significance of the rSZ effect in RX J1347.5-1145 (Zemcov et al.,
2012a).

Recent attempts to measure relativistic SZ effect corrections have been based
on stacking hundreds of galaxy clusters from the Planck all sky survey
(Hurier, 2016; Erler et al., 2018), with these analyses measuring the rSZ
effect signal with a statistical significance similar to what was achieved in
the earlier single-object studies. Looking forward, planned and potential sub-
mm/mm facilities, like LMT/TolTEC (Hughes et al., 2020; Wilson et al., 2020),
the upcoming Fred Young Sub-mm Telescope (FYST; formerly CCAT-prime) (CCAT-
Prime collaboration et al., 2021), the potential Atacama Large Aperture Sub-
mm/mm Telescope (AtLAST; Klaassen et al. 2020a, b) and the Chajnantor
Sub/millimeter Survey Telescope (CSST; Padin 2014) hold the promise of
delivering the required sensitivity, angular resolution, and field of view to
make high-significance rSZ effect measurements routine.

Here we report results from a measurement of the rSZ effect signal toward the
galaxy cluster RX J1347.5-1145 using data from Herschel-SPIRE, Bolocam,
Planck, Chandra and the Hubble Space Telescope (HST). RX J1347.5-1145 is a
famous system that has been the subject of numerous SZ effect studies over the
past two decades (Pointecouteau et al., 1999; Komatsu et al., 2001;
Pointecouteau et al., 2001; Kitayama et al., 2004; Zemcov et al., 2012b;
Sayers et al., 2016a; Kitayama et al., 2016). It is one of the most massive
galaxy clusters observed, with a well-defined cool core and a highly relaxed
morphology over most of its projected surface (Allen et al., 2002; Gitti &
Schindler, 2004; Mantz et al., 2015). However, there is clear evidence for a
shock to the SE of the core, with a corresponding enhancement to the thermal
SZ effect signal in that region (Mason et al., 2010; Plagge et al., 2013; Ueda
et al., 2018; Di Mascolo et al., 2019). Detailed dynamical analyses of RX
J1347.5-1145 indicate that this shock is the result of a 10-to-1 mass ratio
merger occurring primarily in a direction orthogonal to the line of sight
(Johnson et al., 2012). While this merger has produced sloshing features in
the core gas, it has not significantly disrupted the cool core nor the overall
relaxed morphology of the system. Furthermore, the SZ effect enhancement is
relatively small in both amplitude and angular extent compared to the bulk
signal and, as shown in Sayers et al. (2016a), the enhancement is not evident
in the arcminute-resolution imaging we employ in this analysis.

This work has resulted in the development of several novel techniques to
address the challenges of separating unwanted astrophysical contaminants from
the rSZ effect signal. The reduction of the raw data products is described in
Section 2, and our map fitting procedure is outlined in Section 3. A mock
observation pipeline is required to accurately characterize our results, and
the details of this pipeline are described in Section 4. Our analysis of the
mocks to assess the uncertainties and biases in our measured SZ effect
brightness values is reported in Section 5. We present the procedure used to
fit for the average Comptonization parameter and ICM temperature in Section 6.
Finally, we discuss the implications of this work for current and future
studies in Section 7. Throughout this work, we assume a flat $\Lambda$CDM
cosmology with $\Omega_{\rm m}=0.3$, $\Omega_{\Lambda}=0.7,$ and $h=0.7$.

##  2 Data and Reduction

Figure 2: All images show the same 6x6 sq. arcmin field of view. a) SPIRE
three-band false color image. The green box indicates where we fit for the SZ
effect amplitude, is 4′ $\times$ 4′ in size, and is centered on the X-ray
center. b) Chandra X-ray surface brightness image with logarithmic scaling and
smoothing applied. c) HST CLASH optical image with lensing model in red/blue
contours. d) Bolocam SZ effect image. The circle indicates $R_{2500}$, within
which we compute the average SZ effect signal brightness. e) SPIRE three-band
false color image showing a contour overlay of the Bolocam image. The red
cross denotes the X-ray center. Given the astrometric uncertainties, the peak
of the SZ effect in the Bolocam image is coincident with the X-ray center.

###  2.1 SPIRE

SPIRE was an imaging photometer consisting of three focal plane arrays with
approximate band centers at 600, 850, and 1200 GHz, which will be identified
as PLW, PMW and PSW respectively for the remainder of this manuscript Griffin
et al. (2010). The point spread functions (PSFs) are well approximated as
Gaussians, with full-width at half maxima (FWHMs) of 36″, 25″, and 18″. RX
J1347.5-1145 was observed to a depth corresponding to instrument noise levels
of $\lesssim 2$ mJy beam-1 for all three bands as part of the _Herschel_
Multi-Tiered Extragalactic Survey (Oliver et al., 2010, 2012).

The SPIRE science archive data used in this analysis were processed through
the SMAP/SPIRE-HerMES Iterative Map Maker (SHIM) following the description in
Levenson et al. (2010). SHIM is favored over the Herschel-provided tools
because it is optimized to separate large-scale correlated noise from the
signal, making it better suited for the study of extended emission including
the SZ effect. The transfer functions for diffuse astronomical signals were
estimated using the methods described in Viero et al. (2013). The typical
difference between the input and output maps due to high-pass filtering is
$\lesssim 1$ per cent. This is sub-dominant to the $\sim 5$ per cent absolute
calibration uncertainty for SPIRE (see Section 5), so we assume the map-space
signal transfer function is unity (other than the overall mean signal level,
which is not preserved by SHIM).

###  2.2 Bolocam

We use Bolocam photometric imaging data collected at an SZ-emission weighted
band center of 139 GHz with an approximately Gaussian PSF with a FWHM of
59.2″. The square images are 14′ in size, have an overall calibration
uncertainty of 1.7 per cent, and astrometry accurate to 5″. To remove
atmospheric fluctuations from the data, a template subtraction and high-pass
filter are applied to the data timestreams. This results in an effective
angular high-pass filtering of the SZ effect signal, which has been accurately
characterized. To obtain an estimate of the SZ effect surface brightness, we
follow the general method described by Sayers et al. (2019, and references
therein), which we briefly summarize here.

First, an elliptical generalized NFW (gNFW) profile (Nagai et al., 2007),
using the power-law exponents from Arnaud et al. (2010), is fit to the
combination of the Bolocam data and the Planck MILCA $y$-map (Sayers et al.,
2016b; Planck Collaboration et al., 2016a). As noted in Sayers et al. (2016b),
the addition of the Planck data significantly improves the quality of this fit
by constraining the SZ effect signal on large angular scales. The PSF of each
instrument, along with the Bolocam filtering, are fully accounted for in these
fits. Next, the fitted elliptical gNFW model is used to generate a
2-dimensional angular template of the SZ effect signal. This template is then
fit to the Bolocam data, with its normalization as the only free parameter, to
obtain the average surface brightness of the SZ effect signal within R2500 at
a frequency of 139 GHz. R2500 corresponds to the spherical radius enclosing an
average density 2500 times the critical density of the universe, and in RX
J1347.5-1145 is measured to be $0.71\sim$ Mpc based on the analysis of Czakon
et al. (2015). We note that the technique used to obtain R2500 in Czakon et
al. (2015) was based on a generalized scaling relation intended to be
applicable to large, heterogeneous galaxy cluster samples. For highly relaxed
objects, like RX J1347.5-1145, more accurate methods are available (see, e.g.,
Mantz et al., 2014, who find a value of 0.80 Mpc). However, given that R2500
is employed in this work solely as a convenient aperture size that is well-
matched to the observational data, and that Czakon et al. (2015) had
previously computed aperture photometry values from Bolocam data within
R2500$\sim=0.71$ Mpc, we retain that value for this analysis.

In performing the above fits, we assume the map-space noise in the Bolocam
data can be described using a diagonal covariance matrix (i.e., there is no
correlated noise between map pixels). Since this assumption is imperfect, we
characterize the uncertainty on the SZ effect surface brightness by performing
analogous fits to a set of 100 noise realizations. These noise realizations
are identical to those used by Sayers et al. (2019), with the addition of a
kinematic SZ effect signal (see Section 5.2). We find that the distribution of
SZ effect surface brightness values obtained from the fits to these noise
realizations is approximately Gaussian, and so we assign an uncertainty to the
measured SZ effect surface brightness in the observed data based on the
standard deviation of this distribution.

###  2.3 Chandra

ICM temperatures measured using the rSZ effect are approximately pressure-
weighted (e.g., Kay et al., 2008), and so we compute an analogous temperature
from the available Chandra X-ray spectroscopic imaging 111For an example of
the biases that can result from comparing the standard emission-weighted X-ray
spectroscopic temperature to rSZ effect temperatures, see Lee et al. (2020)..
The procedure for reducing and cleaning the data is described in Mantz et al.
(2014, 2015), although we use more recent versions of the Chandra analysis
software and calibration files for this analysis (CIAO version 4.9 and CALDB
version 4.7.4). From these data, we obtain deprojected density and temperature
profiles using the techniques detailed in Mantz et al. (2014, 2016). Using
these profiles, we then compute the pressure-weighted mean temperature of the
ICM within a cylindrical volume defined by R2500 in the plane of the sky and
the maximum radial extent probed by the X-ray data along the line of sight,
corresponding to $\sim 1.8$ R2500. As part of this calculation, we also apply
the empirically-derived temperature calibration bias from Wan et al. (2021),
equal to $0.09\pm 0.13$. With this calibration correction, we obtain a value
of $\langle\textrm{T}_{\textrm{x,pw}}\rangle_{2500}$ $=17.4\pm 2.3$ keV.

###  2.4 HST

The galaxy cluster RX J1347.5-1145 was extensively imaged with _Hubble_ as
part of the Cluster Lensing and Supernova HST (CLASH) survey (Postman et al.,
2012; Koekemoer et al., 2011). The HST images were used to construct mass
models according to the PIEMDeNFW formalism (Zitrin et al., 2015). This
parametric model incorporates elliptical NFW dark matter profiles to model the
cluster dark matter halo, and double pseudo isothermal elliptical mass
distributions to model cluster member galaxies. In order to construct a model
that can be used over the full extent of the SPIRE maps, the best-fit
parameters from the PIEMDeNFW model were used to extrapolate onto a larger
grid (for more details, see Sayers et al., 2019). These grids have sizes of
16′ $\times$ 16′ and a angular resolution of 0.25″.

##  3 Multi-Component Map Fitting

A significant challenge to measuring the SZ effect from observations in the
sub-mm/mm is constraining a signal that is contaminated with emission from
dusty star-forming galaxies (DSFGs), including background DSFGs that are
lensed by the cluster, galactic cirrus dust emission and potentially other
components such as diffuse dust emission from the cluster itself (Vogelsberger
et al., 2019; Erler et al., 2018; Planck Collaboration et al., 2016b, c).
These spatially coincident signals can be disentangled through prior
information on the spatial morphology and underlying spectral energy
distributions (SEDs) of the emission components.

Spectral information has previously been used through multiband matched
filtering to minimize the contribution from unwanted signals in a stacking
analysis of galaxy clusters in Planck maps (Erler et al., 2018). However, in
general Planck does not resolve the galaxy clusters, and so the possibility of
using spatial information is not available. SPIRE data, however, permit
detailed spatial-spectral modeling to separate the unresolved point-like DSFGs
from the diffuse SZ effect and cirrus dust emission. The combination of
angular resolution, instrument sensitivity and the DSFG luminosity function
set intrinsic limits on the efficacy of this type of approach, but for certain
observations the gains in robustness and constraining power can be
significant.

The multi-component fitting method employed in this work builds on the
hierarchical modeling framework of probabilistic cataloging, or PCAT (Brewer
et al., 2013; Daylan et al., 2017; Portillo et al., 2017; Feder et al., 2020).
PCAT is designed to explore the space of catalog models consistent with an
observed image by fitting a Poisson mixture model directly to the data, where
each component represents a point-like source. As a Bayesian hierarchical
model, probabilistic cataloging combines prior information about the source
population of interest with the data likelihood to estimate a posterior
distribution of point-like source models consistent with our prior
expectations and the data. This is represented as an ensemble of catalogs that
naturally encode the often complicated model uncertainties that can arise in
confusion-limited observations (such as from SPIRE) and crowded fields in
general.

The major utility of PCAT as it applies to this analysis is that a point
source model can be estimated directly from the data in such a way that
marginalization over the point source model enables better estimation of a
correlated signal of interest, where in this case the correlated signal is the
SZ effect. This includes marginalization over uncertainties associated with
point source positions and flux densities, as well as uncertainty due to our
ignorance of the true number of DSFGs (down to a given flux density) in the
observed field. This is necessary in the absence of ancillary data from deeper
observations with finer angular resolution. In this section, we describe the
forward model used in probabilistic cataloging and the extension to
observations of the SZ effect toward galaxy clusters.

###  3.1 Generative model

For an image with dimension $(\textrm{W},\textrm{H})$ observed in band $b$,
the surface brightness sampled by pixel $(i,j)$ is written as a sum over
point-like and diffuse components convolved with the PSF $\mathcal{P}$:

$\begin{split}\lambda_{ij}^{b}=\mathcal{P}^{b}\circledast\Bigl{[}\sum_{n=1}^{\textrm{N}_{\textrm{src}}}S_{n}^{b}\delta(x_{i}-x_{n}^{b},y_{j}-y_{n}^{b})\;+\\\
\textrm{A}_{b}^{\textrm{SZ}}I^{\textrm{SZ}}_{b}(x_{i},y_{j})\>+\>B_{ij}^{b}\Bigr{]}.\end{split}$
(1)

The first sum term denotes the contribution from N${}_{\textrm{src}}$ point-
like sources with flux densities
$\\{S_{n}\\}_{n=1}^{\textrm{N}_{\textrm{src}}}$ and positions
$\\{x_{n}\\}_{n=1}^{\textrm{N}_{\textrm{src}}}$, and assumes that the galaxies
in our image are well represented as point sources, which is reasonable given
the SPIRE PSF size and the galaxy redshift distribution. The second term is
the spatially extended SZ effect, included in our model as an angular template
$I_{SZ}^{b}$ scaled by amplitude $\textrm{A}_{b}^{\textrm{SZ}}$, with the
angular template constructed from the elliptical gNFW fit to Bolocam and
Planck (see Section 2.2). The last term, $B_{ij}^{b}$, captures additional
diffuse components in a non-parametric fashion through the addition of a two-
dimensional truncated Fourier series:

$B_{ij}^{b}=B_{0}^{b}+\sum_{n_{x}=1}^{N_{m}}\sum_{n_{y}=1}^{N_{m}}\boldsymbol{\beta}_{n_{x}n_{y}}\cdot\boldsymbol{\mathcal{F}}_{ij}^{n_{x}n_{y}},$
(2)

where $\boldsymbol{\mathcal{F}}_{ij}^{n_{x}n_{y}}$ is a vector of values from
Fourier components with wavevector
$(k_{x},k_{y})=(\textrm{W}/n_{x},\textrm{H}/n_{y})$ evaluated at pixel
$(i,j)$:

$\boldsymbol{\mathcal{F}}_{ij}^{n_{x}n_{y}}=\left(\begin{array}[]{c}\sin\left(\frac{n_{x}\pi
x_{j}}{\textrm{W}}\right)\sin\left(\frac{n_{y}\pi y_{j}}{\textrm{H}}\right)\\\
\sin\left(\frac{n_{x}\pi x_{j}}{\textrm{W}}\right)\cos\left(\frac{n_{y}\pi
y_{j}}{\textrm{H}}\right)\\\ \cos\left(\frac{n_{x}\pi
x_{j}}{\textrm{W}}\right)\sin\left(\frac{n_{y}\pi y_{j}}{\textrm{H}}\right)\\\
\cos\left(\frac{n_{x}\pi x_{j}}{\textrm{W}}\right)\cos\left(\frac{n_{y}\pi
y_{j}}{\textrm{H}}\right)\end{array}\right)$ (3)

and $\boldsymbol{\beta}_{n_{x}n_{y}}$ are the component amplitudes. Including
all four components in $\boldsymbol{\mathcal{F}}^{n_{x}n_{y}}$ is necessary
because we do not wish to impose boundary conditions on the diffuse model. We
determine an appropriate truncation scale for the Fourier component model
empirically by minimizing scatter on the inferred SZ effect template
amplitudes recovered from mock observations (see Section 4). In general, the
choice of truncation scale depends on both the power spectrum of the diffuse
signal contamination and the scale of the PSF, which acts as a low-pass filter
through its convolution with $B_{ij}^{b}$. The mean surface brightness of each
image, captured by $B_{0}^{b}$, is not physically meaningful because the SPIRE
maps are not absolutely calibrated, and so $B_{0}^{b}$ is treated as a
nuisance parameter in our model.

###  3.2 Data likelihood and priors

Our data likelihood is assumed to be Gaussian and is written in map space as a
product over pixels:

$\mathcal{L}=\prod_{b=1}^{\textrm{B}}\prod_{i=1}^{\textrm{W}}\prod_{j=1}^{\textrm{H}}\frac{1}{\sqrt{2\pi\vec{\sigma}^{2}_{b}}}\exp\left(-\frac{(\vec{d}_{b}-\vec{\lambda}_{b})^{2}}{2\vec{\sigma}^{2}_{b}}\right).$
(4)

where $\vec{d}_{b}$ is the observed data vector and $\vec{\lambda}_{b}$ is the
generated model. The noise model that sets the pixel-wise variance
$\sigma_{b}^{2}$ in Eq. (4) is estimated from the SPIRE timestream data
following the procedure described in Viero et al. (2013). The above expression
then reduces to

$\log\mathcal{L}\approx\sum_{b=1}^{\textrm{B}}\sum_{i=1}^{\textrm{W}}\sum_{j=1}^{\textrm{H}}-\frac{(\vec{d}_{b}-\vec{\lambda}_{b})^{2}}{2\vec{\sigma}_{b}^{2}}.$
(5)

We assume galaxies are randomly distributed on the sky, placing uniform priors
on source positions, although we note this is an imperfect assumption mainly
due to the spatial inhomogeneities resulting from gravitational lensing of CIB
sources by the galaxy cluster. The multiband flux density prior $\pi$ is
factorized into a flux density prior for the shortest wavelength band at
250$\mu$m (PSW) and color priors for the remaining two bands:

$\pi(\vec{S})=\pi(S_{\textrm{PSW}})\pi\left(\frac{S_{\textrm{PMW}}}{S_{\textrm{PSW}}}\right)\pi\left(\frac{S_{\textrm{PLW}}}{S_{\textrm{PSW}}}\right).$
(6)

We assume the source flux density distribution follows a power law:

$\pi(S_{\textrm{PSW}})\propto\left(\frac{S_{\textrm{PSW}}}{S^{\textrm{min}}_{\textrm{PSW}}}\right)^{-\alpha}$
(7)

and $\alpha=3.0$ (Casey et al., 2014). The color priors are modeled as
Gaussian with means consistent with the typical DSFG SED. The widths of the
color priors are optimized along with other hyperparameters using the mock
observations described in Section 4 to minimize scatter in the inferred SZ
effect brightness. The minimum source flux density permitted by the model is
determined in a similar fashion and is fixed to $S_{\textrm{PSW}}^{\rm min}=5$
mJy. The diffuse cirrus model is represented in image space as a linear
combination of templates (one template per Fourier component), and the
coefficients of those templates are sampled with uniform priors. Lastly, while
the SZ effect increment between PMW and PLW has a well-defined range of colors
for plausible temperatures, we choose to place independent, uniform priors on
the SPIRE SZ effect template amplitudes to capture potential systematic
effects that could bias the inferred surface brightness values.

By sampling the space of models consistent with observed data $D$ using Eq.
(4) and regularizing the set of solutions with suitable priors, we can compute
the posterior distribution of astrophysical models, $P(\mathcal{M}|D)$,
through Bayes’ rule:

$P(\mathcal{M}|D)=\frac{\mathcal{L}(D|\mathcal{M})\pi(\mathcal{M})}{P(D)}\propto\mathcal{L}(D|\mathcal{M})\pi(\mathcal{M}).$
(8)

Here, $\mathcal{L}(D|\mathcal{M})$ is the likelihood of observing data $D$
given astrophysical signal model $\mathcal{M}$. PCAT uses a Markov Chain Monte
Carlo (MCMC) sampler which has been optimized to efficiently explore the
posterior distribution of catalogs consistent with image data. Details on
implementation and the sampling algorithm can be found in Portillo et al.
(2017) and Feder et al. (2020). The extension of probabilistic cataloging that
performs joint reconstruction with diffuse signal components, PCAT-DE, uses
the same Metropolis-Hastings algorithm to sample mean signal level and
template components in addition to a point source model.

###  3.3 Fitting procedure

Our procedure to fit the model to multiband SPIRE image data happens in two
steps. First, PCAT-DE is run on the PSW-only data to determine a spatial model
for cirrus dust emission. A sixth-order Fourier component model is fit to the
data ($\theta_{\textrm{min}}\sim 1.6\arcmin$, $n_{param}=144$) along with the
point source model and mean surface brightness level. Second, PCAT-DE is run
on all three bands simultaneously, with the shape of the Fourier component
model fixed. The Fourier component model from PSW is then scaled to PMW and
PLW, assuming a constant cirrus SED across the field of view and constant
$B_{0}^{b}$. The assumed cirrus SED is taken from SPIRE observations of the
H-ATLAS SDP field (Bracco et al., 2011). In this step the SZ templates for PMW
and PLW are added to and fit jointly with the rest of the model.

This two-step procedure is performed in order to mitigate the effect of
degeneracies between the cirrus model and the other signal components. The
spectrum of cirrus dust is well constrained and is brightest at short
wavelengths, so determining the spatial structure of cirrus using PSW data
alone is sufficient. The single-band fit also avoids possible bias of the
cirrus dust model by the SZ effect itself, which is also spatially extended.
The SZ effect template amplitudes are floated for PMW and PLW, but since the
the SZ effect brightness is close to zero at high frequencies we fix
$\textrm{A}_{\textrm{PSW}}^{\textrm{SZ}}$ to the value predicated from the
combination of the measured Bolocam SZ effect brightness $\langle
d\textrm{I}_{\textrm{Bcam}}\rangle_{2500}$ and Chandra
$\langle\textrm{T}_{\textrm{x,pw}}\rangle_{2500}$, with $\langle
d\textrm{I}_{\textrm{PSW}}\rangle_{2500}$ =0.005 MJy/sr. This eliminates the
possibility of CIB and/or cirrus emission in the PSW map from being
incorrectly modeled as SZ effect signal.

Probabilistic cataloging is computationally expensive compared to other point
source detection/extraction algorithms. This is because both the number of
parameters in PCAT’s forward model and the degrees of freedom in the image
data are large. To ensure we recover a well-sampled posterior on $\\{$
$\textrm{A}_{b}^{\textrm{SZ}}$ $\\}$, we restrict our multi-component map
fitting procedure to the $4\arcmin\times 4\arcmin$ region centered on the SZ-
defined cluster centroid. However, we compute the best fit mean level for each
band using the larger $10\arcmin\times 10\arcmin$ maps and then fix those
values in fits of the $4\arcmin\times 4\arcmin$ maps. We confirm through tests
on mock data that this procedure does not bias our SZ estimates, but we do
include an additional statistical uncertainty determined by how much the
inferred SZ template amplitudes vary using a range of mean surface brightness
levels consistent with our $10\arcmin\times 10\arcmin$ map fits (see §5).
Crucially, fits on the larger maps are done using the full model from Eq. (1),
which allows for unbiased recovery of the mean level in each band.

The outputs of this fitting process comprise: a catalog of three-band point
source flux densities and associated errors, three-band Fourier component
amplitudes of the cirrus emission, and SZ effect template amplitudes in PMW
and PLW bands expressed as a fitted brightness amplitude
$\textrm{A}_{b}^{\textrm{SZ}}$. These surface brightnesses can then be
corrected for the relevant biases and combined with the measured Bolocam
surface brightness to constrain the SZ effect parameters $\langle
y\rangle_{2500}$ and $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$, as
detailed in Section 6.2.

##  4 Mock Observation Pipeline

In order to validate our analysis pipeline, assess biases on the inferred SZ
effect signal, and quantify uncertainties associated with instrument noise and
astrophysical contaminants, synthetic multiband SPIRE maps of RX J1347.5-1145
are generated. The components in these maps include instrument noise, the SZ
effect, diffuse dust emission from galactic cirrus, and random realizations of
the CIB that include the effects of gravitational lensing. In the following
subsections we describe how these various components are generated. The
individual and combined signal components for the three SPIRE bands are shown
in Figure 3 for one mock observation.

Figure 3: An example of one mock SPIRE observation of the galaxy cluster RX
J1347.5-1145 in PSW (left), PMW (middle) and PLW (right).

###  4.1 Instrument Noise

We use the model described in §2.1 to generate random instrument noise
realizations for our mock observations. Specifically, we assume the noise
fluctuations to be Gaussian and uncorrelated between map pixels, with per
pixel uncertainties determined by the detectors’ time stream variance and
integration time in each pixel. The spatial scan pattern of SPIRE results in a
non-uniform integration time across the field of view that is largest in the
central region and decreases towards the edges (Oliver et al., 2012).

###  4.2 Thermal SZ effect

To include the SZ effect signal in our mock images, we use the elliptical gNFW
model described in §2.2 as a spatial template. The normalization of this
template in the SPIRE bands is computed using SZpack based on the Bolocam
measured surface brightness at 139 GHz and the Chandra-measured pressure-
weighted temperature. The resulting images are then convolved with the
appropriate PSF for each SPIRE band.

###  4.3 Diffuse foregrounds

Initial fits to observed data in RX J1347.5-1145 with PCAT, which assumed a
mean signal level and not a more general diffuse component, produced image
residuals with spatial structure on large angular scales. Data from IRAS
(Miville-Deschênes & Lagache, 2005) and Planck (Planck Collaboration et al.,
2016d, e) show emission with a consistent spatial and spectral signature,
implying that the residual seen in our images is due primarily to thermal dust
emission from galactic cirrus. However, the Planck beam size is very coarse
compared to our map size, meaning only the largest modes are captured. Indeed,
the _Planck_ -interpolated maps were not sufficiently accurate as spatial
templates to model the diffuse emission in the SPIRE maps. As PCAT has a
relatively low minimum source flux density threshold, residual diffuse
emission can be misattributed to low-significance point sources. This was also
observed using _Planck_ -interpolated templates, motivating the Fourier
component approach described in Section 3.1.

Random realizations of cirrus emission are created by drawing Gaussian random
fields with angular power spectra consistent with the observed cirrus signal
in RX J1347.5-1145. In particular, the power spectrum is assumed to scale as
$P(k_{\theta})\propto k_{\theta}^{-2.6}$ (Bracco et al., 2011), with a
normalization set by the measured amplitude of the large-scale power spectrum
of the Planck dust template from Planck Collaboration et al. (2016e),
interpolated to the SPIRE footprint and extrapolated to SPIRE frequencies
using the Planck-estimated parameters of a modified blackbody SED. As with the
other astronomical signals included in our mock images, the resulting
realizations are convolved with the PSF appropriate to each SPIRE band.

###  4.4 CIB

The CIB, which is due primarily to DSFGs but also active galactic nuclei
(AGN), is the brightest astrophysical source of emission at high galactic
latitudes (Hauser & Dwek, 2001). The depth of these observations and the
angular resolution of SPIRE cause the image-space pixel fluctuations to be
dominated by faint, undetected CIB sources rather than the instrument noise,
commonly referred to as “confusion noise” (Nguyen et al., 2010). Further, the
vast majority of bright galaxies in the SPIRE images are not associated with
the cluster, and are instead located behind it (Rawle et al., 2014, 2016).
Since RX J1347.5-1145 is an efficient gravitational lens (Bradac, 2004; Zitrin
et al., 2015), most of the CIB sources in the SPIRE image have been deflected
and magnified. As a result, it is necessary to consider not just the bright
sources that can be detected individually, but also the undetected CIB
sources, many of which have been lensed by the galaxy cluster.

We create mock observations of the CIB in the SPIRE bands using a two-step
process wherein: bright sources individually detected by PCAT are used to
produce constrained CIB realizations and random realizations of the population
of faint undetected sources are generated from an empirical model of the CIB,
including the effects of gravitational lensing on sources behind RX
J1347.5-1145. We describe these steps in detail below.

First, at flux densities above $2\sigma_{\textrm{conf}}$, corresponding to
11.6, 12.6, and 13.6 mJy in PSW, PMW, and PLW (Nguyen et al., 2010), PCAT
detects sources with a completeness of $\gtrsim 90$% and a false detection
rate of $\lesssim 10$% (Feder et al. 2021, in prep.). We thus expect sources
above these thresholds in the PCAT catalog to accurately describe the observed
sky. Therefore, to create a single mock observation of the CIB, we extract the
positions and flux densities of all catalog sources brighter than these
thresholds in at least one SPIRE band. By populating the mocks with different
realizations from PCAT’s catalog ensemble, we ensure that the ensemble of
mocks encodes the measured uncertainties from blind source
extraction/photometry of the CIB.

Second, we generate random catalogs of the positions, SPIRE flux densities,
and redshifts of a set of CIB galaxies using the empirical model from
Béthermin et al. (2012). All of the sources behind the galaxy cluster are
gravitationally lensed according to the mass model derived in Section 2.4. To
retain the correct CIB population statistics, we then remove all of the
sources brighter than $2\sigma_{\textrm{conf}}$ in at least one SPIRE band
from the lensed source catalog. For a single mock, a full catalog of CIB
sources is thus obtained from the combination of bright sources in the PCAT
catalog and fainter sources remaining in the Béthermin et al. (2012) catalog
after this removal. The SPIRE PSF is then used to generate a mock observation
from this catalog.

We note that the empirical model of Béthermin et al. (2012) was updated in
Béthermin et al. (2017, hereafter B12 and B17 respectively), and we also
created mock observations of the CIB using the more recent model. However,
when the effects of gravitational lensing were included, we found a
statistically significant excess of bright sources in mock observations
created from the B17 model compared to the observed data in the PSW and PMW
bands (see Figure 4). Note that the B12 and B17 catalogs were lensed in an
identical manner, and that the catalog generation process is independent of
the lensing step.

The underlying cause of this excess appears to be the gravitational lensing of
numerous faint background sources by the galaxy cluster. However, it is
unclear whether the excess is specific to the geometry of this particular
gravitational lens, or if it generically appears towards massive galaxy
clusters. Regardless of the specific cause of this excess when using the B17
model, it is not observed in the mock observations created from the B12 model,
which show good agreement with the observed data. Therefore, we have used the
older model for this analysis.

Figure 4: Pixel histograms of observed SPIRE maps and the mock observations
created from both the B12 and B17 models. While the mock observations from B12
(orange) show agreement with SPIRE observations, the B17 model (green)
produces an excess of bright pixels.

##  5 Sources of Measurement Error

In this section, the components of the mock observations described in Section
4 are used to assess the impact of each type of signal (or noise) on the
derived SZ results. We also isolate and study the effect of lensing on our
inference by testing CIB realizations both with and without application of the
lensing model.

###  5.1 Instrument noise

For Bolocam, we use the formalism described in Sayers et al. (2011) to produce
random realizations that include both detector noise and fluctuations in the
atmospheric emission. We note that the detector noise has an approximately
flat spectrum, while the atmospheric fluctuations have a power spectrum that
increases as a power law at low angular frequencies in the map.

For SPIRE, the contribution from instrument noise is isolated by generating
random realizations from the SPIRE noise model, adding a fiducial SZ effect
signal to each realization and fitting the signal template to the data within
our analysis pipeline. This may be interpreted as an estimate of the raw
sensitivity of the measurement, in the absence of other systematics.
Quantitative estimates of the uncertainties on the $\langle
d\textrm{I}_{b}\rangle_{2500}$ due to instrument noise are given in Table 1.

###  5.2 Astrophysical contamination - Bolocam

The Bolocam images contain a small, but non-negligible, amount of signal from
unwanted astrophysical contaminants. As originally described in Sayers et al.
(2011), our noise model includes random realizations of the primary CMB
anisotropies and CIB based on their measured angular power spectra. In
addition, the signal from the bright AGN in the BCG is modeled and removed
according to the procedures detailed in Sayers et al. (2013). For this
analysis, we also add a contribution to the noise model due to the kinematic
SZ effect signal resulting from the (unknown) bulk line of sight velocity of
the galaxy cluster. Following the convention of Mueller et al. (2015), we
assume a random velocity centered on zero with a standard deviation of 300 km
sec-1 based on the simulations of Sheth & Diaferio (2001). In computing the
signal from this velocity, we assume the ICM is isothermal with a temperature
equal to the X-ray measured value of 17.4 keV. Quantitative estimates of the
overall uncertainties on $\langle d\textrm{I}_{\textrm{Bcam}}\rangle_{2500}$
arising from these signals are given in Table 1.

###  5.3 Astrophysical contamination - SPIRE

The individual galaxies comprising the CIB contaminate our measurement of the
SZ effect, particularly the spatially correlated emission arising due to
gravitational lensing of the background population. Using the 100 mock CIB
realizations described in Section 4.4, we estimate the associated uncertainty
on the $\langle d\textrm{I}_{b}\rangle_{2500}$ from the aggregate posterior
for the values of $\textrm{A}_{b}^{\textrm{SZ}}$. Because the ensemble of CIB
realizations is well approximated as a collection of independent, identically
distributed draws from an underlying luminosity function, the aggregate
posterior from these mocks should capture the effect of instrument noise, per-
realization CIB model uncertainties, and any error due to intrinsic scatter
from cosmic variance in the CIB. While the B12 model does not include a
clustering term, the fluctuations due to Poisson noise from the CIB dominate
at the scales of interest to our map, and so is assumed to be negligible
(Viero et al., 2013, see Figure 9).

To quantify the contribution of cirrus dust contamination to our error budget,
we compute the scatter on the inferred $\langle d\textrm{I}_{b}\rangle_{2500}$
from an ensemble of 100 lensed mocks without cirrus, and then compare against
the same mock observations with cirrus included. This constrains the effect of
diffuse cirrus emission to the upper limits given in Table 1.

Unlike in the Bolocam data, emission from the AGN in the BCG is relatively dim
compared to the other signals in the SPIRE bands. Specifically, its brightness
is approximately equal to the instrument noise per beam and almost two orders
of magnitude dimmer than the integrated SZ effect signal within R2500 based on
the SED fits from Komatsu et al. (1999) and Sayers et al. (2013). Therefore,
we have not attempted to specifically model the AGN emission in the SPIRE
data.

####  5.3.1 Gravitational lensing

Gravitational lensing has a significant impact on our analysis of the CIB and
the resulting SZ effect constraints. Uncertainties on the $\langle
d\textrm{I}_{b}\rangle_{2500}$ obtained from performing our analysis on lensed
CIB mock realizations are presented in Table 1, and equal to 0.014 and 0.012
MJy sr-1 for PLW and PMW. In addition, we measure a significant bias in the
value of the $\langle d\textrm{I}_{b}\rangle_{2500}$, equal to $-0.019$ and
$-0.021$ MJy sr-1 in the two bands. For comparison, we also performed our
analysis using unlensed CIB mock realizations, obtaining slightly higher
$\langle d\textrm{I}_{b}\rangle_{2500}$ uncertainties of 0.016 and 0.014 MJy
sr-1 and significantly smaller (and positive) biases equal to $+0.007$ MJy
sr-1 in both bands.

These differences between the lensed and unlensed mocks are due primarily to
the effect of “depletion” noted by Blain (2002), where a lack of CIB emission
is observed within the strong lensing region near the center of the galaxy
cluster (Zemcov et al., 2013; Sayers et al., 2019). The effective subtraction
of bright individual sources from the images using PCAT further enhances this
depletion. Thus, the level of CIB emission coincident with the SZ effect
signal is slightly lower than what is observed in unlensed regions of the sky.
In addition, the inferred $\langle d\textrm{I}_{b}\rangle_{2500}$ are biased
low due to the on-average deficit of CIB emission coincident with the SZ
effect signal. For our final analysis, we have corrected the values of
$\langle d\textrm{I}_{b}\rangle_{2500}$ for the measured bias, which, as noted
above, is due primarily to lensing.

###  5.4 Instrumental flux calibration

The SPIRE data are first calibrated according to the procedure described by
Bendo et al. (2013). This calibration is then adjusted using the empirical
cross-calibration factors between Planck-HFI and SPIRE determined by
Bertincourt et al. (2016). As detailed by Bertincourt et al. (2016), the
Planck-HFI calibration has a statistical uncertainty of 1.1 and 1.4 per cent
at 545 and 857 GHz, along with an absolute uncertainty of $\simeq 2$ per cent
and $\simeq 5$ per cent, where the former is obtained from a measurement of
the first two acoustic peaks in the primary CMB anisotropy power spectrum and
the latter is obtained from the ESA2 planetary model of Uranus and the ESA3
planetary model of Neptune. In translating the Planck-HFI calibration to
SPIRE, there is an additional $\simeq 4$ per cent uncertainty due to the SPIRE
PSF calibration, along with a sub-per cent statistical uncertainty in the
cross calibration. Adding these terms in quadrature, we estimate the absolute
SPIRE calibration to be accurate to 4.6 per cent and 6.6 per cent for PLW and
PMW. Furthermore, based on this calibration scheme, we expect a negligible
correlation in the calibration uncertainty between PLW and PMW.

The Bolocam data are calibrated according to the procedure described in Sayers
et al. (2019), which is accurate to 1.7 per cent. In brief, the empirical
model derived in Sayers et al. (2012) is used to correct for variations in
atomospheric transmission based on observing conditions. The planetary model
of Griffin & Orton (1993), corrected based on the Planck-HFI planetary
brightness measurements of Planck Collaboration et al. (2017), is then
employed to determine the absolute calibration.

###  5.5 Additional instrumental calibrations

While there are additional sources of potential systematic errors related to
instrumental calibration (e.g., the measured spectral bandpasses, PSF shape
measurements, astrometric corrections), all of these are likely to be sub-
dominant to the flux calibration accuracy. In addition, many of these
potential sources of systematic error have already been subsumed into the flux
calibration model, and are thus largely accounted for. We therefore do not
explicitly include them in our overall error budget.

Table 1: Sources of Measurement Error

Source |  $<$ Uncertainty on $\langle d\textrm{I}_{b}\rangle_{2500}$ ${}^{\textrm{a}}$ | Notes
---|---|---
| (MJy sr-1) |
Bolocam
Instrument+Atmosphere |  $<$ 0.008${}^{\textrm{a}}$ | See Sayers et al. (2011)
CMB+CIB |  $<$ 0.003${}^{\textrm{a}}$ | See Sayers et al. (2011)
Kinematic SZE |  $<$ 0.007${}^{\textrm{a}}$ | Assuming $\pm 300$ km sec-1 (Sheth & Diaferio, 2001)
Absolute Calibration |  $<$ 0.003${}^{\textrm{a}}$ |  $\pm 1.7$ per cent, see Sayers et al. (2019)
Total |  $<$ 0.012${}^{\textrm{a}}$ | Measured $\langle d\textrm{I}_{\textrm{Bcam}}\rangle_{2500}$ $=-0.153$ MJy sr-1
SPIRE PLW
Instrument |  $<$ 0.003${}^{\textrm{a}}$ | Mock SZ + Instrument Noise run through pipeline
CIB |  $<$ 0.014${}^{\textrm{a}}$ | Mean bias of $-0.019$ MJy sr-1 (primarily due to lensing${}^{\textrm{b}}$)
Cirrus |  $<$0.003 | Mocks with and without cirrus have consistent posteriors
Absolute Calibration |  $<$ 0.003${}^{\textrm{a}}$ |  $\pm 4.6$ per cent based on Bertincourt et al. (2016)
Model fitting |  $<$ 0.002${}^{\textrm{a}}$ | Fixed mean signal level
Total |  $<$ 0.014${}^{\textrm{a}}$ | Measured $\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ $=0.104$ MJy sr-1
SPIRE PMW
Instrument |  $<$ 0.004${}^{\textrm{a}}$ | Mock SZ + Instrument Noise run through pipeline
CIB |  $<$ 0.012${}^{\textrm{a}}$ | Mean bias of $-0.021$ MJy sr-1 (primarily due to lensing${}^{\textrm{b}}$)
Cirrus |  $<$0.003 | Mocks with and without cirrus have consistent posteriors
Absolute Calibration |  $<$ 0.002${}^{\textrm{a}}$ |  $\pm 6.6$ per cent based on Bertincourt et al. (2016)
Model fitting |  $<$ 0.003${}^{\textrm{a}}$ | Fixed mean signal level
Total |  $<$ 0.013${}^{\textrm{a}}$ | Measured $\langle d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ $=0.037$ MJy sr-1

Note. — Estimated uncertainty on the value of $\langle
d\textrm{I}_{b}\rangle_{2500}$ due to both instrumental and astrophysical
sources (see text for more details). The uncertainties due to CIB fluctuations
in the SPIRE PLW and PMW bands are highly correlated, see Figure 5.

${}^{\textrm{a}}$While the CIB-related uncertainty on $\langle
d\textrm{I}_{b}\rangle_{2500}$ is dominated by cosmic variance, a small amount
is also related to fluctuations in the value of N${}_{\textrm{src}}$ in PCAT
(corresponding to 0.004 MJy sr-1 for PLW and 0.003 MJy sr-1 for PMW).

${}^{\textrm{b}}$If the effects of gravitational lensing are not included in
the mocks, then the uncertainty on $\langle d\textrm{I}_{b}\rangle_{2500}$
attributed to variations in the CIB is slightly larger (0.016 MJy sr-1 for PLW
and 0.014 MJy sr-1 for PMW). In addition, the magnitude of the mean bias on
$\langle d\textrm{I}_{b}\rangle_{2500}$ is significantly reduced and its value
is slightly positive ($+0.007$ MJy sr-1 for both PLW and PMW). Thus, the on-
average decrease in CIB brightness near the galaxy cluster center due to
gravitational lensing strongly biases the inferred $\langle
d\textrm{I}_{b}\rangle_{2500}$ to lower values while also minimally decreasing
the associated uncertainty.

##  6 Results

###  6.1 PCAT SPIRE results

The quantities $\langle d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ and $\langle
d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ are inferred using a combination of
the observed SPIRE data and the set of constrained mock observations. To
ensure full coverage over the SZ posterior, we run 100 randomly initialized
Markov chains on the observed data in parallel. Each chain is run for $5\times
10^{6}$ samples, which are then thinned by a factor of 1000. The PCAT-DE model
provides a reasonably good but imperfect description of the data, with typical
reduced-$\chi^{2}$ values, computed using the pixel-wise log-likelihood in Eq.
(5) and the number of parameters at each step of the chain, in the range of
1.2–1.3.

Running several independent chains allows us to assess the level of
convergence in our model. We check that our chains have completed the “burn-
in” phase by visually inspecting the ensemble of trace plots and confirming
that they are well mixed. We discard the first 50% of samples and combine the
remaining ones from all chains to produce the aggregate posterior distribution
of $\textrm{A}_{b}^{\textrm{SZ}}$ values. The mean auto-correlation lengths
from the collection of observed chains are
$\langle\tau\rangle_{\textrm{PMW}}=6$ (thinned) samples and
$\langle\tau\rangle_{\textrm{PLW}}=22$, which results in effective sample
sizes large enough that Monte Carlo error on the results are $<1\%$. We
further validate that our chains are well mixed using the Gelman-Rubin
statistic, also known as the potential scale reduction factor (PSRF) (Gelman,
A., Rubin, D., 1992). For the 100 randomly initialized chains run on the
observed data, $\hat{R}=1.14$ for PMW and $\hat{R}=1.2$ for PLW, which
suggests sufficient convergence to use the aggregated samples from all
observed chains to compute an unbiased estimate of the posterior variance. By
running several chains on individual mock observations, we measure
$\hat{R}<1.1$ and so we do not include any additional correction to the 100
chains run on the mocks, where only one chain is assigned to each mock
realization.

Figure 5 shows the marginalized posterior on
$\textrm{A}_{\textrm{PMW}}^{\textrm{SZ}}$ and
$\textrm{A}_{\textrm{PLW}}^{\textrm{SZ}}$, along with the inferred number of
point-like sources. Both the mock and observed posteriors on
$\textrm{A}_{b}^{\textrm{SZ}}$ in Fig. 5 are uncorrected for the bias due
primarily to lensing depletion (see Section 5.3.1). In both the observed data
and mock realizations, $\textrm{A}_{\textrm{PMW}}^{\textrm{SZ}}$ and
$\textrm{A}_{\textrm{PLW}}^{\textrm{SZ}}$ are positively correlated, with
Pearson correlation coefficients of +0.18 and +0.52 respectively. Since the SZ
effect is independently observed across different SPIRE bands, the correlation
between the two is due to the presence of correlations in the unresolved
source population that are degenerate with the SZ effect signal. In the
observed data, $\textrm{A}_{\textrm{PMW}}^{\textrm{SZ}}$ and
N${}_{\textrm{src}}$ appear uncorrelated. However,
$\textrm{A}_{\textrm{PLW}}^{\textrm{SZ}}$ is anti-correlated with the number
of detected CIB sources with a correlation coefficient of $-0.25$. This
comports with the coarser PLW angular resolution, which results in the CIB
being more spatially degenerate with extended emission from the SZ effect. The
anti-correlation between $\textrm{A}_{\textrm{PLW}}^{\textrm{SZ}}$ and
N${}_{\textrm{src}}$ is washed out in the aggregate mock posterior, where
intrinsic scatter from cosmic variance in the CIB is dominant.

Probing the cross-model covariance between the SZ effect and the union of
point source models with varying N${}_{\textrm{src}}$ is straightforward in
the framework of probabilistic cataloging, since for each SZ posterior sample
there is an associated catalog of CIB sources. To assess the impact of our
imperfect knowledge of N${}_{\textrm{src}}$, we take the quadrature difference
of the fully marginalized uncertainties on $\textrm{A}_{b}^{\textrm{SZ}}$ with
those from the conditional uncertainty assuming the median inferred
N${}_{\textrm{src}}$, i.e.
$\sqrt{\sigma^{2}(\textrm{A}_{b}^{SZ})-\sigma^{2}(\textrm{A}_{b}^{SZ}|\textrm{med}(\textrm{N}_{\textrm{src}}))}$.
By this metric, our results suggest that imperfect knowledge of
N${}_{\textrm{src}}$ results in an uncertainty on $\langle
d\textrm{I}_{b}\rangle_{2500}$ that is $\sim 4$ times smaller than the
uncertainty due to cosmic variance in the CIB.

Additional uncertainties due to the use of a fixed mean level estimated from
the larger $10\arcmin\times 10\arcmin$ maps, inflating the constraints by and
the convergence of the Markov chains are added as a random Gaussian component
to the set of corrected samples, inflating the constraints on $\langle
d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ and $\langle
d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ by 17% and 20%, respectively. After
correcting for the lensing depletion bias, we measure the SZ effect brightness
to be $\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ $=0.104\pm 0.014$
MJy/sr and $\langle d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ $=0.037\pm
0.013$ MJy/sr, corresponding to 7.4$\sigma$ and 2.9$\sigma$ detections of the
SZ effect, respectively.

Finally, the difference between the posterior SZ effect samples and the
_Chandra_ +Bolocam SZ effect amplitude injected into the mocks defines a bias
distribution, which is used to propagate estimates from the observed (biased)
value of $\textrm{A}_{b}^{\textrm{SZ}}$ to the underlying $\langle
d\textrm{I}_{b}\rangle_{2500}$ (see Figure 6 and also the discussion in
Section 5.3.1). To test for model dependence in the bias, mock realizations
with different injected SZ effect amplitudes were analyzed. For amplitudes
spanning the range $[0,2]$ times the nominal _Chandra_ +Bolocam value, the
mean bias is found to be constant within our measurement precision. Thus, it
is valid to apply the linear bias correction to any measured values of
$\textrm{A}_{b}^{\textrm{SZ}}$ within this range.

Figure 5: Marginalized posterior from 100 PCAT chains run on SPIRE maps of RX
J1347.5-1145 for SZ effect template amplitudes
$\textrm{A}_{\textrm{PLW}}^{\textrm{SZ}}$ and
$\textrm{A}_{\textrm{PMW}}^{\textrm{SZ}}$, and the inferred number of CIB
sources N${}_{\textrm{src}}$. Template amplitudes represent the average
surface brightnesses within $R_{2500}$ and have units of MJy sr-1. The SZ
effect posteriors presented here are uncorrected for lensing depletion bias.
The red crosshairs indicate the level of injected SZ effect signal used as
“ground truth” in our mock simulations to estimate the lensing depletion bias,
which is taken from the measured _Chandra_ temperature of
$\langle\textrm{T}_{\textrm{x,pw}}\rangle_{2500}$ $=17.4$ keV combined with
the measured Bolocam $\langle d\textrm{I}_{\textrm{Bcam}}\rangle_{2500}$.
Solid contours cover the 68 and 95 per cent credible regions of each posterior
distribution. Figure 6: 68 and 95 per cent credible regions for
$\textrm{A}_{\textrm{PMW}}^{\textrm{SZ}}$ and
$\textrm{A}_{\textrm{PLW}}^{\textrm{SZ}}$. The posteriors are composed from
the distribution of template amplitudes sampled by PCAT-DE. Credible regions
are shown for the observed RX J1347.5-1145 SPIRE data (purple) and mock SPIRE
maps with instrument noise only (black), instrument noise + unlensed CIB
(green) and instrument noise + lensed CIB (blue). The red star indicates the
injected SZ effect surface brightness values in the mock observations.

###  6.2 SZ spectral likelihood analysis

To obtain constraints on the average Comptonization parameter $\langle
y\rangle_{2500}$ and ICM temperature
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$, we perform a likelihood
analysis of our observed $\langle d\textrm{I}_{b}\rangle_{2500}$ over a range
of $\langle y\rangle_{2500}$ and
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ values, with
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ varying between $[0,60]$ keV.
The model $\langle d\textrm{I}_{b}\rangle_{2500}$ values for each $\langle
y\rangle_{2500}$/$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ grid pair are
computed using SZpack (Chluba et al., 2012, 2013). In all cases, the effective
band center of each $\langle d\textrm{I}_{b}\rangle_{2500}$ used in this
calculation was obtained by weighting the measured spectral bandpass of the
instrument by the temperature-dependent shape of the SZ effect spectrum at the
given value of $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$. The likelihood
assumes independence between the Bolocam and SPIRE measurements of $\langle
d\textrm{I}_{b}\rangle_{2500}$, which is a good approximation given that the
CIB and cirrus signals are negligible in the Bolocam image and the primary CMB
anisotropies and kinematic SZ effect signal are negligible in the SPIRE
images.

Our results using PCAT-DE demonstrate that the inferred distribution of
$\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ and $\langle
d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ values is non-Gaussian to such a
degree that approximating the PDF with a Gaussian covariance matrix
significantly distorts the final constraints. Instead, when computing the
posterior $p(y,T|d\textrm{I})$ with Bayes’ rule:

$p(y,T|d\textrm{I})\propto p(d\textrm{I}|y,T)p(y,T),$ (9)

the likelihood $p(d\textrm{I}|y,T)$ is evaluated directly from the gridded,
bias-corrected samples of $\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$
and $\langle d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ (also referred to as
the “empirical PDF”). The prior $p(y,T)$ places constraints on the range of
spectra that are consistent with our underlying model. In our current
analysis, no prior is placed on the amplitudes of the SZ templates, nor on
their color, at the map-fitting level. While this helps in identifying and
correcting for systematic effects, it does mean a subset of samples from the
posterior fall outside the range of $\langle d\textrm{I}_{b}\rangle_{2500}$
computed by SZPack for $0\leq$ $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$
$\leq 60$ keV. These samples, which reside on the tails of the PDF, comprise
only $\sim 7\%$ of the full sample and are given zero weight in the final
constraints.

The measured $\langle d\textrm{I}_{b}\rangle_{2500}$ and constrained set of SZ
effect spectra are shown in Figure 7, and the posterior on $\langle
y\rangle_{2500}$ and $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ is shown
in Figure 8. Looking separately at the contributions from Bolocam and SPIRE,
one can see that the individual measurements suggest different constraints on
$\langle y\rangle_{2500}$ and $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$.
The Bolocam data constrain the parameters along a positively correlated axis
in the $\langle
y\rangle_{2500}$/$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ plane. The
$\langle y\rangle_{2500}$/$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$
constraints using the empirical SZ posterior from SPIRE data have a more
complicated structure. The preference of a smaller value of $\langle
y\rangle_{2500}$ compared to Bolocam is largely driven by the value of
$\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$, which is approximately
$\sim 2\sigma$ lower than the value expected from the combination of Bolocam
and Chandra. While the Bolocam data are in good agreement with the maximum _a
posteriori_ model estimate, the delta log-likelihood from the SPIRE points is
$\Delta\ln\mathcal{L}=-2.4$, indicating tension with the spectral model.

The joint posterior contains the most probability mass near the tails of the
SPIRE-only posterior, and when combined with Bolocam data leads to bimodal
constraints. This is due to the positive covariance between $\langle
d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ and $\langle
d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ combined with the measured values
falling on opposite sides of the SZ effect spectrum for that range of $\langle
y\rangle_{2500}$/$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$. That is, the
measured values scatter in an anti-correlated manner relative to the model
while the covariance matrix implies a positive correlation, which leads to two
local peaks in $p(d\textrm{I}|y,T)$. There is mild preference for the mode
with larger values of the pair.

Figure 7: Top: SZ effect spectrum for the Bolocam+SPIRE MAP $\langle
y\rangle_{2500}$ and $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ values
(thin blue), along with the spectrum obtained from the Bolocam $\langle
d\textrm{I}_{\textrm{Bcam}}\rangle_{2500}$ and Chandra
$\langle\textrm{T}_{\textrm{x,pw}}\rangle_{2500}$ (red). Also shown is the
best-fit non-relativistic SZ effect spectrum (dashed orange). The blue shaded
region bounds the range of SZ effect spectra within the 68 per cent credible
region of $\langle y\rangle_{2500}$ and
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$. The $\langle
d\textrm{I}_{b}\rangle_{2500}$ values are shown as black points, with reported
errors. Bottom: Residual spectra after subtraction of the non-relativistic SZ
effect spectrum. Note that, due to the negative sign of the SZ effect in the
Bolocam band, a measured $\langle d\textrm{I}_{\textrm{Bcam}}\rangle_{2500}$
below the MAP spectrum indicates a larger rather than a smaller SZ effect
amplitude. The PSW band center is at approximately 1200 GHz, where there is
very little SZ effect signal. Figure 8: The joint posterior distribution
(right) is computed as the product of the SPIRE posterior (left), which is
sampled directly with PCAT-DE, and the Bolocam likelihood (middle) which is
assumed to be Gausisan-distributed with respect to both parameters. The black
circle + yellow cross indicates the position of the maximum _a posteriori_
(MAP) estimate of $\langle y\rangle_{2500}$ and
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$. Log-probabilities are plotted
for visualization purposes. Figure 9: 1D marginalized probability distribution
for cluster temperature $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$, with
joint MAP temperature (black) and credible intervals shaded in blue. Table 2:
Summary of Results

Parameter | MAP value | 68 per cent credible interval | 95 per cent credible interval
---|---|---|---
$\langle y\rangle_{2500}$ $[\times 10^{-4}]$ | 1.56 |  $1.42<$ $\langle y\rangle_{2500}$ $<1.63$ |  $1.29<$ $\langle y\rangle_{2500}$ $<1.71$
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ [keV] | 22.4 |  $10.4<$ $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ $<33.0$ |  $0.0<$ $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ $<39.5$

Note. — Fitted values of $\langle y\rangle_{2500}$ and
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ for RX J1347.5-1145.

Our maximum a posteriori (MAP) estimate for the ICM temperature is
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$$=22.4$ keV, with a 68% highest
posterior density credible interval (HDPI) of $10.4<$
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ $<33.0$ keV (see Figure 9).
This estimate is consistent with the X-ray predicted, pressure-weighted
temperature of $\langle\textrm{T}_{\textrm{x,pw}}\rangle_{2500}$ $=17.4\pm
2.3$ keV, indicating good agreement between measurements.

To understand the sensitivity of our data in the absence of the bimodality
noted above, we also repeat the same analysis using the $\langle
d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ and $\langle
d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ values determined from the
combination of Bolocam and Chandra (i.e., $\langle
d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ $=0.134$ MJy sr-1 and $\langle
d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ $=0.033$ MJy sr-1). This eliminates
the issue of positively-covariant values falling on opposite sides of the SZ
effect spectrum, and removes the bimodality, while fully preserving the noise
as characterized by the empirical PDF. Under these conditions, the posterior
is single-peaked and the constraints marginalized over $\langle
y\rangle_{2500}$ yield a MAP $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$
$=12.3$ keV with a 68 per cent credible interval of $5.8<$
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ $<20.5$ keV. This represents a
reduction in the range of the credible interval by 35 per cent compared to the
observed data.

If we instead repeat the analysis using the observed values of $\langle
d\textrm{I}_{b}\rangle_{2500}$, but with randomized samples to remove the
correlation between $\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ and
$\langle d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ in the empirical PDF, then
we find a 68 per cent credible interval of $5.1<$
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ $<20.2$ keV. This is nearly
identical to the interval obtained using _Chandra_ +Bolocam predicted $\langle
d\textrm{I}_{b}\rangle_{2500}$ values to eliminate the bimodality, and
suggests that the $\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$/$\langle
d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ correlation similarly degrades our
constraints on $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$.

As illustrated in Table 1, CIB fluctuations dominate the uncertainties on
$\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ and $\langle
d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$. To quantitatively assess the impact
of the CIB on our derived $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$
constraints, we also perform our analysis on the SZ + instrument noise only
mock posterior (see Figure 6), finding a MAP estimate of
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$=17.1 keV with a 68 per cent
credible interval of $12.8<$ $\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$
$<19.0$ keV.

##  7 Discussion

We have combined observations from SPIRE, Bolocam, and Planck to constrain the
thermodynamic properties of the ICM in the galaxy cluster RX J1347.5-1145
through a measurement of the rSZ effect. In order to probe the desired SZ
effect signal, we accounted for significant contamination from unwanted
astrophysical components, in particular the CIB. Not only is the CIB brighter
than the SZ effect signal in the SPIRE bands, but the individual sources that
comprise it have a range of SEDs. In addition, most CIB sources reside behind
RX J1347.5-1145 and have thus been deflected and magnified due to
gravitational lensing. To properly characterize the effects of the CIB and
other relevant signals such as galactic cirrus emission, we developed a mock
observation pipeline to fully assess the uncertainties and biases in measuring
the SZ effect brightness $\langle d\textrm{I}_{b}\rangle_{2500}$. Our mock
pipeline includes a mass model derived from HST, allowing us to accurately
capture the effects of gravitational lensing, which we found to be the largest
source of bias on the derived values of $\langle
d\textrm{I}_{b}\rangle_{2500}$.

Rather than treating each astrophysical source in isolation, we employ an
extension of the forward modeling framework, PCAT-DE, to simultaneously
constrain the signal components. This enables us to make more efficient use of
the available information in the observed data while robustly capturing
correlated, non-Gaussian uncertainties on $\langle
d\textrm{I}_{b}\rangle_{2500}$. We propagate samples from the PCAT-DE
posterior to estimates of $\langle d\textrm{I}_{b}\rangle_{2500}$ and use the
resulting empirical PDF to constrain the Comptonization parameters. This was
important in our analysis – we discovered that approximating the SPIRE
posterior with a Gaussian covariance matrix significantly biased the final
constraints.

As noted in Section 6.2, our constraints on
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ are strongly influenced by the
relatively low measured value of $\langle
d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$. A large positive correlation
between the measured $\langle d\textrm{I}_{\textrm{PMW}}\rangle_{2500}$ and
$\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ values, coupled with them
falling on opposite sides of the SZ effect spectrum for intermediate
temperatures, results in a bimodal posterior for
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ and a significant expansion of
the 68 per cent credible interval range (22.6 keV compared to 14.7 keV when a
more likely value of $\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ is
assumed). We have not been able to identify any significant biases in the
measured $\langle d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$, but it is
possible that one or more do exist. For instance, the particular spatial
distribution of astrophysical signals toward RX J1347.5-1145 may result in
PCAT incorrectly associating some of the PLW SZ effect signal with CIB
sources. It is also possible that some unmodeled astrophysical components are
significant enough to produce a bias, such as the excess SZ effect signal
associated with an ICM shock SE of the cluster center (Komatsu et al., 2001;
Mason et al., 2010; Kitayama et al., 2016; Di Mascolo et al., 2019). RX
J1347.5-1145 also has a pronounced cool core, giving rise to large temperature
variations as a function of radius within the ICM. This invalidates the
assumption of isothermality along the line of sight, and causes a slight
change in the SZ spatial structure between the Bolocam and SPIRE observing
bands. Since we are utilizing Bolocam for the SZ effect spatial template, this
could change the location and significance of the actual SZ effect “peak”.
Another possibility is that the unexplained cluster-coincident emission with a
dust-like SED noted by Erler et al. (2018) has not been fully described by the
combination of the PCAT source catalog and diffuse Fourier component model.

Applying our analysis method to a larger sample of galaxy clusters will allow
us to better assess if the somewhat low measured value of $\langle
d\textrm{I}_{\textrm{PLW}}\rangle_{2500}$ is due to a statistical fluctuation,
something particular to RX J1347.5-1145, or a more general bias in our
methodology. Multi-probe SPIRE, Bolocam, HST, Chandra, and Planck observations
exist for an additional 20 galaxy clusters, and another 18 objects can be
studied using data from the ACT survey in place of Bolocam (Aiola et al.,
2020). In addition to enabling checks for subtle, emergent systematics, this
sample of 39 galaxy clusters will allow us to place tight constraints on the
ensemble-average temperature measured using the rSZ effect. Taking the raw
sensitivity from this study in the absence of a bimodal
$\langle\textrm{T}_{\textrm{sz}}\rangle_{2500}$ posterior, we estimate that
the sample mean temperature will be measured to a precision of approximately 1
keV. By comparing with pressure-weighted X-ray temperatures, an empirical
temperature calibration of Chandra and/or XMM-Newton should thus be possible.
This level of sensitivity would also allow us to probe the expected
differences between rSZ effect and X-ray emission-weighted temperatures due to
spatial and thermal inhomogeneities in the ICM (Lee et al., 2020).

Looking forward, potential future wide-field and high angular resolution
facilities such as AtLAST and CSST will provide better decoupling of the SZ
effect from unwanted astrophysical signals, resulting in higher signal-to-
noise ratio rSZ effect measurements of the ICM temperature while probing
smaller physical scales. These data would complement future X-ray facilities
and open exciting new areas of study, such as routine temperature measurements
in the highest redshift galaxy clusters and resolved imaging of the hottest
shock-heated gas produced in major merger events. In this work, we have
demonstrated an analysis framework that can help maximize the science reach of
these facilities.

## Acknowledgments

The scientific results reported in this article were based on data obtained
from the HerMES project (http://hermes.sussex.ac.uk/), a Herschel Key
Programme utilising Guaranteed Time from the SPIRE instrument team, ESAC
scientists and a mission scientist. The HerMES data was accessed through the
Herschel Database in Marseille (http://hedam.lam.fr) operated by CeSAM and
hosted by the Laboratoire d’Astrophysique de Marseille. The Bolocam data are
publicly available at
https://irsa.ipac.caltech.edu/data/Planck/release_2/ancillary-
data/bolocam/bolocam.html.

The authors acknowledge Research Computing at the Rochester Institute of
Technology for providing computational resources and support that have
contributed to the research results reported in this publication.

The authors would also like to acknowledge Glenn Morris
(https://orcid.org/0000-0003-2985-9962) for their assistance with Chandra data
reduction.

This work was supported by NASA ADAP research grant 80NSSC19K1018. AZ
acknowledges support from the Ministry of Science and Technology, Israel.

## References

  * Aiola et al. (2020) Aiola, S., Calabrese, E., Maurin, L., et al. 2020, J. Cosmology Astropart. Phys, 2020, 047, doi: 10.1088/1475-7516/2020/12/047
  * Allen et al. (2002) Allen, S. W., Schmidt, R. W., & Fabian, A. C. 2002, MNRAS, 335, 256, doi: 10.1046/j.1365-8711.2002.05554.x
  * Arnaud et al. (2010) Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92, doi: 10.1051/0004-6361/200913416
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bendo et al. (2013) Bendo, G. J., Griffin, M. J., Bock, J. J., et al. 2013, MNRAS, 433, 3062, doi: 10.1093/mnras/stt948
  * Bertincourt et al. (2016) Bertincourt, B., Lagache, G., Martin, P. G., et al. 2016, A&A, 588, A107, doi: 10.1051/0004-6361/201527313
  * Béthermin et al. (2012) Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJ, 757, L23, doi: 10.1088/2041-8205/757/2/L23
  * Blain (2002) Blain, A. W. 2002, MNRAS, 330, 219, doi: 10.1046/j.1365-8711.2002.05058.x
  * Bleem et al. (2015) Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27, doi: 10.1088/0067-0049/216/2/27
  * Bracco et al. (2011) Bracco, A., Cooray, A., Veneziani, M., et al. 2011, MNRAS, 412, 1151, doi: 10.1111/j.1365-2966.2010.17971.x
  * Bradac (2004) Bradac, M. 2004, in Baryons in Dark Matter Halos, ed. R. Dettmar, U. Klein, & P. Salucci, 14.  https://arxiv.org/abs/astro-ph/0410642
  * Brewer et al. (2013) Brewer, B. J., Foreman-Mackey, D., & Hogg, D. W. 2013, AJ, 146, 7, doi: 10.1088/0004-6256/146/1/7
  * Béthermin et al. (2017) Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, Astronomy & Astrophysics, 607, A89, doi: 10.1051/0004-6361/201730866
  * Casey et al. (2014) Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45, doi: 10.1016/j.physrep.2014.02.009
  * CCAT-Prime collaboration et al. (2021) CCAT-Prime collaboration, Aravena, M., Austermann, J. E., et al. 2021, arXiv e-prints, arXiv:2107.10364.  https://arxiv.org/abs/2107.10364
  * Chluba et al. (2012) Chluba, J., Nagai, D., Sazonov, S., & Nelson, K. 2012, MNRAS, 426, 510, doi: 10.1111/j.1365-2966.2012.21741.x
  * Chluba et al. (2013) Chluba, J., Switzer, E., Nelson, K., & Nagai, D. 2013, MNRAS, 430, 3054, doi: 10.1093/mnras/stt110
  * Czakon et al. (2015) Czakon, N. G., Sayers, J., Mantz, A., et al. 2015, ApJ, 806, 18, doi: 10.1088/0004-637X/806/1/18
  * Daylan et al. (2017) Daylan, T., Portillo, S. K. N., & Finkbeiner, D. P. 2017, ApJ, 839, 4, doi: 10.3847/1538-4357/aa679e
  * Di Mascolo et al. (2019) Di Mascolo, L., Churazov, E., & Mroczkowski, T. 2019, MNRAS, 487, 4037, doi: 10.1093/mnras/stz1550
  * Erler et al. (2018) Erler, J., Basu, K., Chluba, J., & Bertoldi, F. 2018, MNRAS, 476, 3360, doi: 10.1093/mnras/sty327
  * Fabian et al. (2011) Fabian, A. C., Sanders, J. S., Allen, S. W., et al. 2011, MNRAS, 418, 2154, doi: 10.1111/j.1365-2966.2011.19402.x
  * Feder et al. (2020) Feder, R. M., Portillo, S. K. N., Daylan, T., & Finkbeiner, D. 2020, AJ, 159, 163, doi: 10.3847/1538-3881/ab74cf
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Gelman, A., Rubin, D. (1992) Gelman, A., Rubin, D. 1992, Statistical Science, 7, 4, 457-511
  * Gitti & Schindler (2004) Gitti, M., & Schindler, S. 2004, A&A, 427, L9, doi: 10.1051/0004-6361:200400086
  * Griffin & Orton (1993) Griffin, M. J., & Orton, G. S. 1993, Icarus, 105, 537, doi: 10.1006/icar.1993.1147
  * Griffin et al. (2010) Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3, doi: 10.1051/0004-6361/201014519
  * Hauser & Dwek (2001) Hauser, M. G., & Dwek, E. 2001, ARA&A, 39, 249, doi: 10.1146/annurev.astro.39.1.249
  * Hauser et al. (1998) Hauser, M. G., Arendt, R. G., Kelsall, T., et al. 1998, ApJ, 508, 25, doi: 10.1086/306379
  * Hughes et al. (2020) Hughes, D. H., Schloerb, F. P., Aretxaga, I., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11445, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 1144522, doi: 10.1117/12.2561893
  * Hurier (2016) Hurier, G. 2016, A&A, 596, A61, doi: 10.1051/0004-6361/201629726
  * Itoh et al. (1998) Itoh, N., Kohyama, Y., & Nozawa, S. 1998, ApJ, 502, 7, doi: 10.1086/305876
  * Johnson et al. (2012) Johnson, R. E., Zuhone, J., Jones, C., Forman, W. R., & Markevitch, M. 2012, ApJ, 751, 95, doi: 10.1088/0004-637X/751/2/95
  * Kaiser (1986) Kaiser, N. 1986, MNRAS, 222, 323, doi: 10.1093/mnras/222.2.323
  * Kay et al. (2008) Kay, S. T., Powell, L. C., Liddle, A. R., & Thomas, P. A. 2008, MNRAS, 386, 2110, doi: 10.1111/j.1365-2966.2008.13183.x
  * Kitayama et al. (2004) Kitayama, T., Komatsu, E., Ota, N., et al. 2004, PASJ, 56, 17, doi: 10.1093/pasj/56.1.17
  * Kitayama et al. (2016) Kitayama, T., Ueda, S., Takakuwa, S., et al. 2016, PASJ, 68, 88, doi: 10.1093/pasj/psw082
  * Klaassen et al. (2020a) Klaassen, P. D., Mroczkowski, T. K., Cicone, C., et al. 2020a, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11445, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114452F, doi: 10.1117/12.2561315
  * Klaassen et al. (2020b) Klaassen, P. D., Mroczkowski, T. K., Cicone, C., et al. 2020b, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11445, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114452F, doi: 10.1117/12.2561315
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, The Astrophysical Journal Supplement Series, 197, 36, doi: 10.1088/0067-0049/197/2/36
  * Komatsu et al. (1999) Komatsu, E., Kitayama, T., Suto, Y., et al. 1999, ApJ, 516, L1, doi: 10.1086/311983
  * Komatsu et al. (2001) Komatsu, E., Matsuo, H., Kitayama, T., et al. 2001, PASJ, 53, 57, doi: 10.1093/pasj/53.1.57
  * Lee et al. (2020) Lee, E., Chluba, J., Kay, S. T., & Barnes, D. J. 2020, MNRAS, 493, 3274, doi: 10.1093/mnras/staa450
  * Levenson et al. (2010) Levenson, L., Marsden, G., Zemcov, M., et al. 2010, MNRAS, 409, 83, doi: 10.1111/j.1365-2966.2010.17771.x
  * Mantz et al. (2014) Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2014, MNRAS, 440, 2077, doi: 10.1093/mnras/stu368
  * Mantz et al. (2016) Mantz, A. B., Allen, S. W., Morris, R. G., & Schmidt, R. W. 2016, MNRAS, 456, 4020, doi: 10.1093/mnras/stv2899
  * Mantz et al. (2015) Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2015, MNRAS, 449, 199, doi: 10.1093/mnras/stv219
  * Markevitch et al. (2002) Markevitch, M., Gonzalez, A. H., David, L., et al. 2002, ApJ, 567, L27, doi: 10.1086/339619
  * Markevitch & Vikhlinin (2007) Markevitch, M., & Vikhlinin, A. 2007, Phys. Rep., 443, 1, doi: 10.1016/j.physrep.2007.01.001
  * Mason et al. (2010) Mason, B. S., Dicker, S. R., Korngut, P. M., et al. 2010, ApJ, 716, 739, doi: 10.1088/0004-637X/716/1/739
  * McDonald et al. (2017) McDonald, M., Allen, S. W., Bayliss, M., et al. 2017, ApJ, 843, 28, doi: 10.3847/1538-4357/aa7740
  * McNamara & Nulsen (2007) McNamara, B. R., & Nulsen, P. E. J. 2007, ARA&A, 45, 117, doi: 10.1146/annurev.astro.45.051806.110625
  * Miville-Deschênes & Lagache (2005) Miville-Deschênes, M.-A., & Lagache, G. 2005, ApJS, 157, 302, doi: 10.1086/427938
  * Mueller et al. (2015) Mueller, E.-M., de Bernardis, F., Bean, R., & Niemack, M. D. 2015, ApJ, 808, 47, doi: 10.1088/0004-637X/808/1/47
  * Nagai et al. (2007) Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1, doi: 10.1086/521328
  * Nguyen et al. (2010) Nguyen, H. T., Schulz, B., Levenson, L., et al. 2010, Astronomy and Astrophysics, 518, L5, doi: 10.1051/0004-6361/201014680
  * Oliver et al. (2010) Oliver, S. J., Wang, L., Smith, A. J., et al. 2010, A&A, 518, L21, doi: 10.1051/0004-6361/201014697
  * Oliver et al. (2012) Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614, doi: 10.1111/j.1365-2966.2012.20912.x
  * Padin (2014) Padin, S. 2014, Appl. Opt., 53, 4431, doi: 10.1364/AO.53.004431
  * Plagge et al. (2013) Plagge, T. J., Marrone, D. P., Abdulla, Z., et al. 2013, ApJ, 770, 112, doi: 10.1088/0004-637X/770/2/112
  * Planck Collaboration et al. (2013) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013, A&A, 550, A131, doi: 10.1051/0004-6361/201220040
  * Planck Collaboration et al. (2016a) Planck Collaboration, Aghanim, N., Arnaud, M., et al. 2016a, A&A, 594, A22, doi: 10.1051/0004-6361/201525826
  * Planck Collaboration et al. (2016b) Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016b, A&A, 596, A104, doi: 10.1051/0004-6361/201628522
  * Planck Collaboration et al. (2016c) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016c, A&A, 594, A23, doi: 10.1051/0004-6361/201527418
  * Planck Collaboration et al. (2016d) Planck Collaboration, Aghanim, N., Ashdown, M., et al. 2016d, A&A, 596, A109, doi: 10.1051/0004-6361/201629022
  * Planck Collaboration et al. (2016e) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016e, A&A, 586, A132, doi: 10.1051/0004-6361/201424945
  * Planck Collaboration et al. (2017) Planck Collaboration, Akrami, Y., Ashdown, M., et al. 2017, A&A, 607, A122, doi: 10.1051/0004-6361/201630311
  * Pointecouteau et al. (1999) Pointecouteau, E., Giard, M., Benoit, A., et al. 1999, ApJ, 519, L115, doi: 10.1086/312117
  * Pointecouteau et al. (2001) —. 2001, ApJ, 552, 42, doi: 10.1086/320447
  * Portillo et al. (2017) Portillo, S. K. N., Lee, B. C. G., Daylan, T., & Finkbeiner, D. P. 2017, AJ, 154, 132, doi: 10.3847/1538-3881/aa8565
  * Postman et al. (2012) Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25, doi: 10.1088/0067-0049/199/2/25
  * Pratt et al. (2019) Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25, doi: 10.1007/s11214-019-0591-0
  * Prokhorov & Colafrancesco (2012) Prokhorov, D. A., & Colafrancesco, S. 2012, MNRAS, 424, L49, doi: 10.1111/j.1745-3933.2012.01284.x
  * Prokhorov et al. (2011) Prokhorov, D. A., Colafrancesco, S., Akahori, T., et al. 2011, MNRAS, 416, 302, doi: 10.1111/j.1365-2966.2011.19037.x
  * Rawle et al. (2014) Rawle, T. D., Altieri, B., Egami, E., et al. 2014, MNRAS, 442, 196, doi: 10.1093/mnras/stu868
  * Rawle et al. (2016) —. 2016, MNRAS, 459, 1626, doi: 10.1093/mnras/stw712
  * Reiprich et al. (2013) Reiprich, T. H., Basu, K., Ettori, S., et al. 2013, Space Sci. Rev., 177, 195, doi: 10.1007/s11214-013-9983-8
  * Sarazin (1986) Sarazin, C. L. 1986, Reviews of Modern Physics, 58, 1, doi: 10.1103/RevModPhys.58.1
  * Sayers et al. (2012) Sayers, J., Czakon, N. G., & Golwala, S. R. 2012, ApJ, 744, 169, doi: 10.1088/0004-637X/744/2/169
  * Sayers et al. (2011) Sayers, J., Golwala, S. R., Ameglio, S., & Pierpaoli, E. 2011, ApJ, 728, 39, doi: 10.1088/0004-637X/728/1/39
  * Sayers et al. (2013) Sayers, J., Mroczkowski, T., Czakon, N. G., et al. 2013, ApJ, 764, 152, doi: 10.1088/0004-637X/764/2/152
  * Sayers et al. (2016a) Sayers, J., Zemcov, M., Glenn, J., et al. 2016a, ApJ, 820, 101, doi: 10.3847/0004-637X/820/2/101
  * Sayers et al. (2016b) Sayers, J., Golwala, S. R., Mantz, A. B., et al. 2016b, ApJ, 832, 26, doi: 10.3847/0004-637X/832/1/26
  * Sayers et al. (2019) Sayers, J., Montaña, A., Mroczkowski, T., et al. 2019, ApJ, 880, 45, doi: 10.3847/1538-4357/ab29ef
  * Sheth & Diaferio (2001) Sheth, R. K., & Diaferio, A. 2001, MNRAS, 322, 901, doi: 10.1046/j.1365-8711.2001.04202.x
  * Sunyaev & Zeldovich (1972) Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments on Astrophysics and Space Physics, 4, 173
  * Ueda et al. (2018) Ueda, S., Kitayama, T., Oguri, M., et al. 2018, ApJ, 866, 48, doi: 10.3847/1538-4357/aadd9d
  * Viero et al. (2013) Viero, M. P., Wang, L., Zemcov, M., et al. 2013, ApJ, 772, 77, doi: 10.1088/0004-637X/772/1/77
  * Vogelsberger et al. (2019) Vogelsberger, M., McKinnon, R., O’Neil, S., et al. 2019, MNRAS, 487, 4870, doi: 10.1093/mnras/stz1644
  * Wan et al. (2021) Wan, J. T., Mantz, A. B., Sayers, J., et al. 2021, arXiv e-prints, arXiv:2101.09389.  https://arxiv.org/abs/2101.09389
  * Wilson et al. (2020) Wilson, G. W., Abi-Saad, S., Ade, P., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11453, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 1145302, doi: 10.1117/12.2562331
  * Wright (1979) Wright, E. L. 1979, ApJ, 232, 348, doi: 10.1086/157294
  * Zemcov et al. (2010) Zemcov, M., Rex, M., Rawle, T. D., et al. 2010, A&A, 518, L16, doi: 10.1051/0004-6361/201014685
  * Zemcov et al. (2012a) Zemcov, M., Aguirre, J., Bock, J., et al. 2012a, ApJ, 749, 114, doi: 10.1088/0004-637X/749/2/114
  * Zemcov et al. (2012b) —. 2012b, ApJ, 749, 114, doi: 10.1088/0004-637X/749/2/114
  * Zemcov et al. (2013) Zemcov, M., Blain, A., Cooray, A., et al. 2013, ApJ, 769, L31, doi: 10.1088/2041-8205/769/2/L31
  * Zhao (2003) Zhao, H. 2003, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0312165
  * Zitrin et al. (2015) Zitrin, A., Fabris, A., Merten, J., et al. 2015, The Astrophysical Journal, 801, 44, doi: 10.1088/0004-637x/801/1/44