Imaging the Thermal and Kinematic Sunyaev-Zel'dovich Effect Signals in a Sample of Ten Massive Galaxy Clusters: Constraints on Internal Velocity Structures and Bulk Velocities

TL;DR

This study used multi-wavelength observations to image the Sunyaev-Zel'dovich effect in ten galaxy clusters, constraining their internal velocities and bulk motions, revealing significant internal velocity structures in merging clusters.

Contribution

First to combine SZ, X-ray, and lensing data to map internal velocities and bulk motions in a sample of galaxy clusters with high spatial resolution.

Findings

Detected internal velocity dispersions >1000 km/s in merging clusters.

Measured ensemble-mean line-of-sight velocity of 430 km/s with uncertainties.

Identified correlation between merger orientation and velocity variance.

Abstract

We have imaged the Sunyaev-Zel'dovich (SZ) effect signals at 140 and 270 GHz towards ten galaxy clusters with Bolocam and AzTEC/ASTE. We also used Planck data to constrain the signal at large angular scales, Herschel-SPIRE images to subtract the brightest galaxies that comprise the cosmic infrared background (CIB), Chandra imaging to map the electron temperature $T_e$ of the intra-cluster medium (ICM), and HST imaging to derive models of each galaxy cluster's mass density. The galaxy clusters gravitationally lens the background CIB, which produced an on-average reduction in brightness towards the galaxy clusters' centers after the brightest galaxies were subtracted. We corrected for this deficit, which was between 5-25% of the 270 GHz SZ effect signal within $R_{2500}$. Using the SZ effect measurements, along with the X-ray constraint on $T_e$, we measured each galaxy cluster's average line of sight (LOS) velocity $v_z$ within $R_{2500}$, with a median per-cluster uncertainty of +-700 km/s. We found an ensemble-mean <$v_z$> of 430+-210 km/s, and an intrinsic cluster-to-cluster scatter $\sigma_{int}$ of 470+-340 km/s. We also obtained maps of $v_z$ over each galaxy cluster's face with an angular resolution of 70". All four galaxy clusters previously identified as having a merger oriented along the LOS showed an excess variance in these maps at a significance of 2-4$\sigma$, indicating an internal $v_z$ rms of $\gtrsim$1000 km/s. None of the six galaxy clusters previously identified as relaxed or plane of sky mergers showed any such excess variance.

Tables

Table 1. Table 1: Galaxy Cluster Sample

Name	RA	Dec	Redshift	M500	140 GHz rms	270 GHz rms	Dynamical State
	HH:MM:SS.s	DD:MM:SS		${10}^{14}$ M☉	MJy sr-1	MJy sr-1
Abell 0697	08:42:57.6	$+$36:21:57	0.282	$17.1\pm 2.9$	0.010	0.025 (B)	LOS-merger
Abell 1835	14:01:01.9	$+$02:52:40	0.253	$12.3\pm 1.4$	0.011	0.031 (B)	relaxed
MACS J0018.5$+$1626	00:18:33.4	$+$16:26:13	0.546	$16.5\pm 2.5$	0.013	0.019 (A)	LOS-merger
MACS J0025.4$-$1222	00:25:29.9	$-$12:22:45	0.584	0$7.6\pm 0.9$	0.011	0.025 (A)	POS-merger
MACS J0454.1$-$0300	04:54:11.4	$-$03:00:51	0.538	$11.5\pm 1.5$	0.010	0.024 (A)	POS-merger
MACS J0717.5$+$3745	07:17:32.1	$+$37:45:21	0.546	$24.9\pm 2.7$	0.020	0.020 (B)	LOS-merger
MACS J2129.4$-$0741	21:29:25.7	$-$07:41:31	0.589	$10.6\pm 1.4$	0.015	0.023 (A)	LOS-merger
RX J0152.7$-$1357	01:52:41.1	$-$13:58:07	0.833	0$7.8\pm 3.0$	0.014	0.014 (A)	POS-merger
RX J1226.9$+$3332	12:26:57.9	$+$33:32:49	0.888	0$7.8\pm 1.1$	0.015	0.021 (B)	POS-merger
RX J1347.5$-$1145	13:47:30.8	$-$11:45:09	0.451	$21.7\pm 3.0$	0.013	0.032 (A)	POS-merger

Table 2. Table 2: Instrument Band Centers

Observing Band	Blackbody	Thermal SZ	Kinematic SZ	Synchrotron	Thermal Dust
Bolocam 140 GHz	140.5 GHz	139.3 GHz	140.1 GHz	139.2 GHz	141.2 GHz
Bolocam 270 GHz	270.9 GHz	274.9 GHz	268.0 GHz	267.7 GHz	272.2 GHz
AzTEC 270 GHz	271.3 GHz	275.5 GHz	268.1 GHz	267.8 GHz	272.8 GHz

Table 3. Table 3: SZ Effect Calibration Uncertainty

Observing Band	Measurement	Planck Abs	Extrapolation	PSF	Atm Trans	Total
Bolocam 140 GHz	1.0%	0.6%	0.0%	1.2%	0.2%	1.7%
Bolocam 270 GHz	1.0%	0.7%	1.3%	2.6%	0.3%	3.2%
AzTEC 270 GHz	1.2%	0.7%	1.3%	2.8%	0.3%	3.4%

Table 4. Table 4: Chandra X-ray Observations

Name	Inst.	ObsIDs	Usable Exp. Times (ksec)
Abell 0697	ACIS-I	4217	19.2
Abell 1835	all ACIS-I	6880,6881,7370	115.9,36.3,39.5
MACS J0018.5$+$1626	ACIS-I	520	64.1
MACS J0025.4$-$1222	all ACIS-I	3251,5010,10413,10786,10797	18.0,23.8,75.6,13.7,23.8
MACS J0454.1$-$0300	ACIS-I,ACIS-S	529,902	13.7,41.9
MACS J0717.5$+$3745	all ACIS-I	1655*,4200,16235,16305	—,54.9,67.3,89.9
MACS J2129.4$-$0741	all ACIS-I	3199*,3595	—,18.2
RX J0152.7$-$1357	ACIS-I	913	34.7
RX J1226.9$+$3332	all ACIS-I	3180,5014	29.1,30.8
RX J1347.5$-$1145	all ACIS-I	3592,13516,13999,14407	56.6,39.0,54.4,63.0

Table 5. Table 5: DSFG Detections

Name	270 GHz Det.	270 GHz Lim.	Herschel–SPIRE Det.	600 GHz Lim.	Counterparts
Abell 0697	04	4.00 mJy	121	04.64 mJy	03
Abell 1835	02	4.84 mJy	057	09.16 mJy	01
MACS J0018.5$+$1626	23	2.60 mJy	111	04.42 mJy	19
MACS J0025.4$-$1222	11	3.20 mJy	123	04.68 mJy	09
MACS J0454.1$-$0300	19	2.96 mJy	038	08.62 mJy	09
MACS J0717.5$+$3745	13	3.76 mJy	110	05.58 mJy	13
MACS J2129.4$-$0741	20	3.24 mJy	012	11.98 mJy	08
RX J0152.7$-$1357	22	2.20 mJy	060	06.62 mJy	16
RX J1226.9$+$3332	05	4.52 mJy	027	010.00 mJy	04
RX J1347.5$-$1145	09	3.60 mJy	088	04.52 mJy	09

Table 6. Table 6: The Impact of CIB Lensing on the Measured 270 GHz SZ Effect Brightness

Name	No Lensing Correction	With Lensing Correction	Difference
Abell 0697	0.037 MJy sr-1	0.046 MJy sr-1	0.009 MJy sr-1
Abell 1835	0.069 MJy sr-1	0.080 MJy sr-1	0.011 MJy sr-1
MACS J0018.5$+$1626	0.097 MJy sr-1	0.117 MJy sr-1	0.020 MJy sr-1
MACS J0025.4$-$1222	0.043 MJy sr-1	0.052 MJy sr-1	0.009 MJy sr-1
MACS J0454.1$-$0300	0.081 MJy sr-1	0.091 MJy sr-1	0.010 MJy sr-1
MACS J0717.5$+$3745	0.109 MJy sr-1	0.131 MJy sr-1	0.022 MJy sr-1
MACS J2129.4$-$0741	0.040 MJy sr-1	0.053 MJy sr-1	0.013 MJy sr-1
RX J0152.7$-$1357	0.081 MJy sr-1	0.086 MJy sr-1	0.005 MJy sr-1
RX J1226.9$+$3332	0.105 MJy sr-1	0.124 MJy sr-1	0.019 MJy sr-1
RX J1347.5$-$1145	0.078 MJy sr-1	0.087 MJy sr-1	0.009 MJy sr-1

Table 7. Table 7: Derived ICM Parameters

Name	Temperature (keV)	Optical Depth (${10}^{-3}$)	Bulk $v_z$ (km s-1)	Internal $v_z$ rms (km s-1)
Abell 0697	0${8.99}_{-0.42}^{+0.53}$	0${4.88}_{-0.99}^{+0.86}$	$+{1620}_{-1500}^{+1250}$	${1820}_{-940}^{+940}$
Abell 1835	0${7.66}_{-0.13}^{+0.13}$	0${5.69}_{-1.05}^{+1.04}$	00$-{70}_{-960}^{+850}$	$\le 1970$ (95% CL)
MACS J0018.5$+$1626	0${8.30}_{-0.40}^{+0.49}$	0${7.64}_{-1.01}^{+0.99}$	0$-{110}_{-640}^{+610}$	${810}_{-490}^{+600}$
MACS J0025.4$-$1222	0${5.67}_{-0.24}^{+0.25}$	0${5.83}_{-1.04}^{+1.05}$	0$+{530}_{-750}^{+700}$	$\le 1570$ (95% CL)
MACS J0454.1$-$0300	0${8.83}_{-0.50}^{+0.56}$	0${6.64}_{-1.10}^{+1.10}$	0$+{370}_{-1000}^{+910}$	$\le 1590$ (95% CL)
MACS J0717.5$+$3745	${12.83}_{-1.42}^{+1.42}$	0${7.55}_{-1.00}^{+0.92}$	0$+{740}_{-560}^{+530}$	${1260}_{-360}^{+430}$
MACS J2129.4$-$0741	0${8.52}_{-1.14}^{+1.44}$	0${6.13}_{-1.23}^{+1.16}$	$+{1570}_{-870}^{+810}$	${1170}_{-510}^{+560}$
RX J0152.7$-$1357	0${4.72}_{-0.59}^{+0.56}$	0${9.96}_{-2.24}^{+2.12}$	0$+{150}_{-610}^{+560}$	$\le 960$ (95% CL)
RX J1226.9$+$3332	0${6.84}_{-0.59}^{+0.75}$	${10.13}_{-1.36}^{+1.47}$	00$+{40}_{-480}^{+450}$	$\le 1260$ (95% CL)
RX J1347.5$-$1145	0${9.47}_{-0.29}^{+0.37}$	0${6.94}_{-0.71}^{+0.70}$	0$+{950}_{-680}^{+640}$	$\le 860$ (95% CL)

Equations

$$F(\nu )=F_0(1-e^{(-\nu /{\nu }_0{)}^{\beta }})B(\nu ,T_{\text{d}})$$

Full text

galaxies: clusters: intracluster medium — cosmology: observations

††journal: ApJ

1 Introduction

Velocity measurements have long been used to probe the detailed properties of large-scale structure, for example the velocity dispersion of cluster galaxies (Zwicky, 1937) and the galaxy rotation curves that provided evidence of dark matter (Rubin et al., 1980). As another more recent example, the Hitomi X-ray satellite provided the first direct measurement of the velocity structure of the intra-cluster medium (ICM) in the core of the Perseus cluster (Hitomi Collaboration et al., 2016, 2018), providing new insights on the interaction between the ICM and the central active galactic nucleus as well as a large-scale velocity shear due to cosmic accretion and mergers (Lau et al., 2017; ZuHone et al., 2018). In addition, the statistical properties of the cosmological velocity field can be used to constrain a range of parameters, particularly those related to dark energy and possible modifications of general relativity (e.g., Kaiser 1987; Percival & White 2009). To date, nearly all velocity measurements have been obtained via spectroscopy, mainly at optical wavelengths (e.g., Abolfathi et al. 2018). One challenge to these spectroscopic measurements is the fundamental degeneracy between the object’s recessional velocity due to the expansion of the universe and its peculiar velocity relative to that expansion. The kinematic Sunyaev-Zel’dovich (SZ) effect signal, which is a Doppler shift of cosmic microwave background (CMB) photons inverse Compton scattering with a distribution of electrons, has long held the promise of addressing this challenge by providing velocity measurements relative to the fixed reference frame of the CMB (Sunyaev & Zeldovich 1972, 1980; for a recent review see Mroczkowski et al. 2019).

However, measurements of the kinematic SZ effect signal have proven difficult, mainly due to a lack of raw sensitivity but also due to contamination from a range of unwanted astronomical signals (e.g., Benson et al. 2003, 2004; Kitayama et al. 2004; Zemcov et al. 2012). This situation is slowly changing, as a range of modern instruments have been able to obtain tentative detections of the kinematic SZ effect in resolved observations of exceptional individual galaxy clusters with very high velocity sub-components (Sayers et al., 2013a; Adam et al., 2017) and in aggregate for large statistical samples (Planck Collaboration et al., 2016a; Soergel et al., 2016; De Bernardis et al., 2017). Looking forward, the next generation of instrumentation aims to advance from these first detections to detailed studies using the kinematic SZ effect (e.g., Morandi et al. 2013; Mittal et al. 2018). While these SZ effect studies are unlikely to reach the velocity sensitivity demonstrated by Hitomi in the central regions of nearby galaxy clusters, they will ideally complement future X-ray observations from facilities such as XRISM, Athena, and Lynx by providing velocity measurements at higher redshifts and/or further from the galaxy cluster’s center.

In this work, we used observations from Bolocam and AzTEC/ASTE, along with ancillary data from Herschel–SPIRE, Chandra, Planck, and the Hubble Space Telescope (HST) to obtain resolved images of the SZ effect signal towards a sample of ten galaxy clusters. This analysis was built upon the previous work of Sayers et al. (2013a), who used a subset of these data to detect the kinematic SZ effect signal towards one of the galaxy clusters in our sample, MACS J0717.5$+$3745. In Section 2, we describe the sample of ten galaxy clusters in detail. The datasets and their associated reduction (including the reconstruction of lens models) are then presented in Sections 3 and 4. Our fits to the SZ effect signals, and the galaxy cluster-averaged bulk velocities obtained from these fits are given in Section 5. We then present resolved images of the SZ effect signals in Section 6. Finally, we provide a summary of our analysis in Section 7.

2 Galaxy Cluster Sample

This study was focused on a sample of ten massive galaxy clusters with available data from Bolocam/AzTEC, Herschel-SPIRE, Chandra, and HST. A brief description of the dynamical state of each galaxy cluster is given below, with a summary in Table 1.

1. Abell 0697:

Girardi et al. (2006), based on Chandra X-ray and galaxy cluster member spectroscopic measurements, suggested that this system is undergoing a complex merger mainly along the line of sight (LOS). This complex merger scenario is further supported by the detailed study of its giant radio halo by Macario et al. (2010). Rossetti et al. (2013) also found indications for a merger mainly along the LOS. 2. Abell 1835:

This galaxy cluster was among the first targets of both Chandra and XMM-Newton, and that imaging revealed a highly relaxed morphology (Peterson et al., 2001; Schmidt et al., 2001). A wide range of subsequent studies have supported the conclusion that this is one of the most relaxed known galaxy clusters (e.g., Mantz et al. 2015). 3. MACS J0018.5$+$1626:

Solovyeva et al. (2007) found this galaxy cluster to be undergoing a merger based on Chandra and XMM-Newton data, and Piffaretti et al. (2003) found evidence for LOS elongation based on a joint X-ray and SZ effect analysis. Mann & Ebeling (2012), in their systematic study of 108 galaxy clusters to search for binary mergers, found this galaxy cluster to have a morphological code of 3 on their scale of 1–4, with 4 being the most likely to be undergoing a major merger. However, the reason it was not classified as a 4 was the relatively small offset between the BCG and the X-ray peak, which would be consistent with a merger primarily along the LOS. 4. MACS J0025.4$-$1222:

This galaxy cluster is a dramatic plane of sky (POS) merger, similar to the Bullet Cluster, and has been studied in detail by several groups (Bradač et al., 2008; Ma et al., 2010; Riseley et al., 2017; Cibirka et al., 2018). Mann & Ebeling (2012) listed this galaxy cluster as a textbook example of a binary merger and gave it a morphological code of 4. 5. MACS J0454.1$-$0300:

Both Donahue et al. (2003) and Jeltema et al. (2005) found the X-ray morphology of this galaxy cluster to be elongated in the E–W direction in the POS, indicating a possible merger along that orientation. Furthermore, Mann & Ebeling (2012) gave this galaxy cluster a morphological code of 3, and found a significant offset between the BCG and the X-ray peak. 6. MACS J0717.5$+$3745:

The detailed analysis of Ma et al. (2009) showed this galaxy cluster to be a complex merger with a significant component along the LOS. In particular, they identified four merging subclusters in the system, and they labeled the largest subcluster, which is located slightly SE of the X-ray center, as “C”. Approximately 1.5′ NW of “C” is subcluster “B”, which appears to be moving with a LOS velocity of $+3000$ km s*-1* relative to “C”. This scenario was further supported by a range of subsequent analyses, including two based on kinematic SZ effect measurements (Mann & Ebeling, 2012; Sayers et al., 2013a; Adam et al., 2017; van Weeren et al., 2017). 7. MACS J2129.4$-$0741:

This galaxy cluster was given a morphological code of 3 by Mann & Ebeling (2012), and was described in that paper as a complex merger that is occurring primarily along the LOS. 8. RX J0152.7$-$1357:

Maughan et al. (2006), based on XMM-Newton data, found that this galaxy cluster is undergoing a merger along two main axes, both oriented in the POS. A consistent merger scenario was found by Molnar et al. (2012) based on the offset between the X-ray and SZ effect signal peaks.

9. RX J1226.9$+$3332:

Maughan et al. (2007) found evidence for merger activity in a joint Chandra and XMM-Newton analysis. The weak lensing analysis of Jee & Tyson (2009) further supported a merger scenario. They found a large POS separation of the clumps, indicating that the merger may be oriented primarily along the POS. More recent SZ effect imaging from Korngut et al. (2011) and Adam et al. (2015) provided additional evidence for a POS merger scenario. 10. RX J1347.5$-$1145:

A range of independent analyses have found evidence for a merger in the core region of this galaxy cluster, oriented along the SW–NE direction and primarily in the POS (Mason et al., 2010; Johnson et al., 2012; Plagge et al., 2013; Kreisch et al., 2016; Ueda et al., 2018).

3 Datasets

3.1 Bolocam 140 GHz

All of the galaxy clusters in our sample were imaged with Bolocam at 140 GHz,111Throughout this work we refer to the SZ effect bands as “140 GHz” and “270 GHz”. The precise band centers for a range of source spectra are given in Table 2. and all of those data have been used in previous analyses (e.g., Sayers et al. 2013a; Czakon et al. 2015) and are publicly available.222https://irsa.ipac.caltech.edu/data/Planck/release_2/ancillary-data/bolocam/ The images have a point-spread function (PSF) with a solid angle that corresponds to a Gaussian with a full-width at half maximum (FWHM) of 59.2″. The data were collected in 10 minute observations using a sinusoidal Lissajous pattern with differing periods in the right ascension and declination directions, resulting in a coverage that drops to half its peak value at a radius of 5–6′. We obtained approximately 100 such individual observations per galaxy cluster. Astrometry, with an rms uncertainty of $\simeq 5\arcsec$, was computed based on frequent observations of nearby bright objects.

Nightly observations of Uranus and Neptune were used to calibrate the detector response, and a single empirical fit as a function of atmospheric opacity, accurate to 1.0%, was computed for all of the nights within a given observing run (typically $\sim 10$ nights, see Sayers et al. 2012). For this work, we used the planetary models from Griffin & Orton (1993) rescaled based on the recent measurements from Planck, which are accurate to 0.6% at 140 GHz (Planck Collaboration et al., 2017). While our empirical fit accounted for changes in band-averaged atmospheric transmission as a function of opacity, it did not account for the slight changes in the spectral shape of the atmospheric transmission, which we estimated to produce a 0.2% rms uncertainty in our calibration (see Sayers et al. 2012).

In addition, in transferring the calibration from the point-like planets to resolved SZ effect surface brightness measurements there was an additional uncertainty due to our characterization of the PSF solid angle, which we estimated to be 1.2% based on the quadrature sum of two separate uncertainties. First, Sayers et al. (2009) measured the per-detector solid angle with an rms of 3.1%, with no evidence for variation from detector to detector. Therefore, averaging over the $\simeq 100$ optical detectors resulted in a 0.3% rms measurement uncertainty. Second, the measured solid angle was based on a source spectrum matching that of Uranus and Neptune, which were used for the PSF calibration measurements. We assumed the PSF was diffraction limited, which means its solid angle was different for sources with different spectral shapes, such as the thermal and kinematic SZ effect signals. To account for this difference we included an additional rms uncertainty of 1.2%, equal to the average difference in diffraction-limited PSF solid angle for the effective band centers of the thermal and kinematic SZ effect signals compared to the effective band centers for Uranus and Neptune.

In total, we estimated our calibration to be accurate to an rms uncertainty of 1.7% (see Table 3).

3.2 Bolocam 270 GHz

Four of the galaxy clusters in our sample were observed with Bolocam at 270 GHz, using the same observing strategy detailed above for the 140 GHz data. The 270 GHz Bolocam images have PSFs with a solid angle that corresponds to a Gaussian with a FWHM of 33.2″. Compared to the 140 GHz data, some of the uncertainties on the calibration were slightly different for the 270 GHz data (see Table 3). Specifically, the absolute Planck measurements were accurate to 0.7% and the PSF solid angle characterization resulted in a 2.6% calibration uncertainty (0.6% due to measurement uncertainty and 2.5% due to the differing effective band centers of the thermal and kinematic SZ effect signals). In addition, unlike at 140 GHz, there was no Planck band centered near our observing band at 270 GHz. As a result, we extrapolated the Planck measurements at 220 and 350 GHz, and we estimated this extrapolation resulted in a 1.3% uncertainty based on the deviations obtained from calibrating the Griffin & Orton (1993) model at one of those frequencies and then comparing its prediction to the measured value at the other frequency. The total calibration uncertainty was determined to have an rms uncertainty of 3.2%.

3.3 AzTEC 270 GHz

Six of the galaxy clusters in our sample were observed with AzTEC at 270 GHz from the ASTE telescope (AzTEC was built as a nearly exact replica of Bolocam, see Wilson et al. 2008). The scan pattern used for these observations was very similar to the Lissajous used in the Bolocam observations, and the resulting coverage was similar. The PSF in the images has a solid angle that corresponds to a Gaussian with a FWHM of 30.4″. The calibration uncertainty was nearly identical to the 270 GHz Bolocam data, although the slightly higher measurement uncertainty resulted in a total calibration uncertainty with an rms of 3.4% (see Table 3).

3.4 Herschel–SPIRE

All of the galaxy clusters in our sample were observed by Herschel–SPIRE as part of either the Herschel Multi-tiered Extragalactic Survey (HerMES, Oliver et al. 2012) or the Herschel Lensing Survey (HLS, Egami et al. 2010). Herschel–SPIRE was a three-band photometric imager operating at 600, 850, and 1200 GHz with PSFs with FWHMs of 18.1″, 25.2″, and 36.6″(Griffin et al., 2010). The absolute calibration uncertainty of the Herschel–SPIRE data was 5.5% for unresolved sources, and was verified by cross-calibrating with Planck (Bertincourt et al., 2016). In all cases, the Herschel–SPIRE coverage was sufficient to produce images in all three bands comparable in size to the Bolocam and AzTEC images.

3.5 Chandra

Each galaxy cluster was observed in one or more Chandra X-ray imaging observations. The observation identification numbers (ObsIDs) and exposure times are listed in Table 4. Additionally, we provide information about whether the observation was taken with the imaging or spectroscopic Advanced CCD Imaging Spectrometer (ACIS-I or ACIS-S, respectively). Since both instruments were used in imaging mode, this only impacted the sensitivity, background, and field of view of the exposure. Since each CCD array subtends $8\arcmin\times 8\arcmin$, observations with either ACIS-I or ACIS-S covered a sufficiently large field of view for this analysis.

3.6 HST

We reconstructed lens models for the ten galaxy clusters of our sample using multiband HST imaging, essential for the identification of multiple-image constraints. Although the lens models were largely based on existing models, for completeness we describe the latest available HST imaging which enabled these models. Eight galaxy clusters from our sample were imaged extensively with both optical and near-infrared broadbands in the framework of large lensing surveys such as Cluster Lensing And Supernova with HST (CLASH; PI: Postman; MACS J0717.5$+$3745, MACS J2129.4$-$0741, RX J1347.5$-$1145, RX J1226.9$+$3332, Reionization Cluster Survey (RELICS; PI: Coe; Abell 0697, MACS J0018.5$+$1626, MACS J0025.4$-$1222, RX J0152.7$-$1357), and the HST Frontier Fields (PIs: Mountain, Lotz; MACS J0717.5$+$3745). For the remaining two galaxy clusters, reduced images were downloaded from the HST Legacy Archive, taken in program ID 11591 for both Abell 1835 and MACS J0454.1$-$0300 (PI: Kneib), and programs IDs 10493 (PI: Gal-Yam), 9722 (PI: Ebeling), 9292 (PI: Ford), and 9836 (PI: Ellis), for MACS J0454.1$-$0300. The typical depth for most galaxy clusters was $\sim 26.5-27$ AB per band, and the typical pixel scale was 0.05″–0.06″ per pixel. Details of the lens modeling are given in Section 4.4.

4 Data Reduction

4.1 Bolocam and AzTEC

The Bolocam data at 140 and 270 GHz, along with the AzTEC data at 270 GHz, were reduced in a uniform manner using the analysis pipeline described in detail in Sayers et al. (2011). For the SZ effect analysis, a template of the atmospheric brightness fluctuations was computed by averaging the signal from all of the detectors at each time sample within a single $\simeq 10$ minute observation. A single correlation coefficient between each detector’s data stream and the template was then computed, and the template was subtracted after rescaling by this correlation coefficient. For Bolocam, the correlation coefficient was computed using only the data within a narrow bandwidth of the two fundamental Lissajous scan frequencies. For AzTEC, where the scan frequencies were constantly modulated, we instead computed the correlation coefficient using all of the data within the bandwidth 0.5–2.0 Hz. After this subtraction, a high-pass filter was applied to the data streams, with a characteristic frequency of 250 mHz for the 140 GHz data and 500 mHz for the 270 GHz data. The template removal and high-pass filter resulted in a non-unity transfer function for astronomical signals, and we computed a single transfer function for the two-dimensional image of each galaxy cluster at each observing frequency according to the procedure described in Sayers et al. (2011).

At 270 GHz, for both AzTEC and Bolocam, we also performed a second data reduction using an adaptive principal component analysis (PCA) in place of the average template subtraction (Laurent et al., 2005; Aguirre et al., 2011). The adaptive PCA method was not as effective as the average template subtraction for recovering the SZ effect signal from the galaxy cluster, but it was better for detecting unresolved objects (Sayers et al., 2013a).

Regardless of the subtraction algorithm, the noise properties of the images were estimated using a set of 1000 random realizations based on the procedure given in Sayers et al. (2016a). First, 1000 jackknife realizations were generated by creating images after randomly selecting half of the individual observations and multiplying their data by $-1$. On average, this procedure removed all of the astronomical signals while preserving the noise properties of the instrument and the atmospheric fluctuations. To each of these 1000 jackknife images, a random realization of the primary CMB fluctuations, the background population of dusty star-forming galaxies (DSFGs) that comprise the cosmic infrared background (CIB), and the population of radio galaxies were added. Each instrument’s PSF and subtraction-dependent signal transfer function was accounted for prior to adding these astronomical source realizations. In order to fully capture any correlations in these unwanted astronomical signals between 140 and 270 GHz, we did not generate separate realizations at the two observing frequencies. Instead, a single realization was scaled to both frequencies.

After producing the images, along with their associated noise realizations, we then jointly fitted an elliptical generalized NFW (gNFW) model (Nagai et al., 2007a) to the 140 GHz Bolocam images and the Planck all-sky $y$–map (Planck Collaboration et al., 2016b), according to the method detailed in Sayers et al. (2016b) which fully accounted for the Bolocam transfer function and the Planck and Bolocam PSFs. For these fits, the normalization and scale radius of the model were varied while fixing the three power law exponents $\alpha$, $\beta$, and $\gamma$ to the best fit values of Arnaud et al. (2010). For the radial scales typically probed by our data, \simeq 0.3$$R_{2500}–3.0$R_{2500}$,333 The values of $R_{2500}$ used in this work were taken from Czakon et al. 2015, where they were computed from Chandra X-ray data using a scaling relation between gas mass and total mass. this model had sufficient freedom to provide a good fit quality (see Sayers et al. 2011 and Czakon et al. 2015), particularly since ellipticity in the POS was allowed. Furthermore, while a range of more recent observational studies have found different best-fit values of $\alpha$, $\beta$, and $\gamma$ (e.g., Sayers et al. 2013a; Planck Collaboration et al. 2013; Ghirardini et al. 2019), the actual profile shapes are in excellent agreement owing to the strong degeneracies between the parameters, particularly when the scale radius is allowed to vary, as it was in our analysis. We therefore do not expect any significant biases due to our choice of model to describe the shape of the SZ effect signal.

The resulting best-fit gNFW model was then subtracted from the adaptive-PCA-reduced 270 GHz images, accounting for the transfer function and PSF of those images. This subtraction removed most of the SZ effect signal, leaving the background CIB as the dominant astronomical signal in the images. We then used StarFinder to detect all of the unresolved objects with a $\textrm{S/N}>4$ from the resulting images. We typically detected $\simeq 10$ such objects in each image, all of which were presumed to be DSFGs (see Table 5). As detailed below in Section 4.2, Herschel–SPIRE was more sensitive to the signal from DSFGs, and typically detected an order of magnitude more objects.

For the next step in our analysis, we returned to the 140 and 270 GHz SZ effect images created using the average template subtraction. From these images, we subtracted all of the radio galaxies listed in Sayers et al. (2013c) and all of the DSFGs detected in the 270 GHz images and/or the Herschel–SPIRE images. To subtract the radio galaxies, the power law fits from Sayers et al. (2013c) were extrapolated to 140 and 270 GHz. The DSFGs were categorized into three groups, with a slightly different procedure used to subtract the sources from within each of these groups. The first group included DSFGs detected at 270 GHz without a counterpart identified in the Herschel–SPIRE detections. These were subtracted from the 270 GHz data based on their detected flux density, and from the 140 GHz data based on a rescaling of the flux density according to $\nu^{2.5}$. The second group included DSFGs detected at 270 GHz which had a Herschel–SPIRE counterpart. For these sources, the 270 GHz and Herschel–SPIRE three-band measurements were simultaneously fitted to a greybody SED of the form:

[TABLE]

where the values of the normalization $F_{0}$ and the dust temperature $T_{\textrm{\footnotesize d}}$ were varied, $\nu$ is the observed frequency, $\nu_{0}=3000$ GHz (Draine, 2006), $\beta=1.95$ is the dust emissivity spectral index,444 Given the noise and spectral coverage of the data, we were unable to robustly constrain the values of both $\beta$ and $T_{\textrm{\tiny d}}$ for a single source. We therefore fixed the value of $\beta$ in our fits. To determine what value of $\beta$ to use, we compared the measured flux density in the 270 GHz images at the position of every Herschel–SPIRE detection to the flux density computed from a greybody fit solely to the three-band Herschel–SPIRE images extrapolated to 270 GHz for a range of fixed $\beta$ values. On average, the two flux densities agree within the measurement noise for $\beta=1.95\pm 0.11$. This was consistent with the value of $\beta$ found in several other recent studies (e.g., Magnelli et al. 2012; Smith et al. 2013). and $B(\nu,T_{\textrm{\footnotesize d}})$ is the Planck function. The sources were then subtracted from the 140 and 270 GHz images based on the flux densities obtained from this greybody fit. The third and final group included DSFGs detected by Herschel–SPIRE that were not associated with a 270 GHz detection. These sources were subtracted in an analogous way to those in the second group, except that the greybody SED was fit solely to the three-band Herschel–SPIRE images.

This procedure for characterizing and subtracting the DSFGs was nearly identical to what was described in detail in the appendix of Sayers et al. (2013a). Since the quality of the data used in this work was nearly identical to the data used by Sayers et al. (2013a), the same overall implications were also true and are summarized here. In particular, all sources brighter than $\simeq 1$ mJy at 270 GHz were detected, and some sources were detected down to a limit of $\simeq 0.1$ mJy at 270 GHz. In aggregate, these detected sources represent $\simeq 30$% of the total emission from the CIB at that frequency. As noted above, most of the detections were made by Herschel–SPIRE, and the AzTEC/Bolocam detections typically amounted to only $\simeq 5$–10% of the total CIB. Even after subtracting $\simeq 30$% of the CIB emission, the fluctuations due to the remaining sources added an rms per beam of approximately 0.5 mJy, and these fluctuations degraded our SZ effect constraints at 270 GHz by $\simeq 10$–20% compared to what would have been possible with perfect removal of the CIB. While these undetected CIB sources added a non-negligible amount of noise, they did not produce a measurable bias in the SZ effect constraints, likely because their distribution was well described by a Gaussian rms given the PSF size and noise level typical of our 270 GHz data.

In addition, we subtracted an image of the average apparent signal deficit in the CIB produced by galaxy cluster lensing of the DSFG population when the brightest individual sources were removed. This effect was first detected by Herschel–SPIRE, and was used to estimate the total brightness of the CIB (Zemcov et al., 2013). In addition, Lindner et al. (2015) measured a lower than expected SZ effect signal in Herschel–SPIRE, and they speculated that this was due to lensing of the CIB based on the previous Zemcov et al. (2013) results. To estimate the lensing-induced CIB deficit in our SZ effect images, we propagated a random realization of the CIB through the lensing model determined from the HST data (see Section 4.4).

For each galaxy cluster we generated 100 such realizations. Individual bright sources were then removed from these realizations in a way that mimicked the procedure applied to the actual data, which resulted in the detection limits given in Table 5. After removing the bright sources, each realization was then spatially filtered based on the transfer function for the data reduction using an average template subtraction. The 100 realizations were then averaged for each galaxy cluster, and the result was subtracted from the actual SZ effect images. While lensing can produce large brightness variations in the CIB due to the (rare) high magnification of intrinsicly bright DSFGs, all such extremely bright objects were subtracted from both our real data and the 100 lensed realizations. As a result, the typical brightness fluctuations between the 100 lensed realizations were well described by the unlensed CIB realizations already included in our noise model. Examples of the average lensed CIB are shown in Figure 1. Based on the bulk SZ effect fits described in Section 5, the typical deficit in the CIB due to lensing was $\simeq 15$% of the SZ effect brightness at 270 GHz (see Table 6).

4.2 Herschel–SPIRE

The three-band Herschel–SPIRE images were used to search for and characterize DSFG candidates. The data were reduced using the Herschel Interactive Processing Environment (HIPE, Ott et al. 2006; Ott 2010) and the HerMES SMAP package (Levenson et al., 2010; Viero et al., 2013). A list of DSFG candidates was compiled based on the SCAT procedure (Smith et al., 2012), with the requirement that each source have a $\textrm{S/N}>2$ at both 600 and 850 GHz. We found that many of the brighter DSFGs at 270 GHz were not detected at 1200 GHz by Herschel–SPIRE, and so we did not impose a S/N threshold on those data. This typically resulted in $\simeq 100$ DSFG candidates per galaxy cluster (see Table 5).

4.3 Chandra

The Chandra data reduction and analysis closely followed the methods presented in Ogrean et al. (2015) and van Weeren et al. (2017), based on the publicly available scripts used in those previous analyses.555See https://github.com/gogrean/MACS-J0717-Filament/blob/master/code/notebooks/. Briefly, the data were reprocessed to apply the latest calibration at the time, in this case CIAO 4.10 with CALDB 4.7.8. Both of these tools were released sufficiently after each observation used for analysis that the calibration was stable/unchanging for newer releases. In the case of observations taken in VFAINT mode, the check_vf_phaevents option was used to provide additional filtering for background events. As in Ogrean et al. (2015), we extracted light curves from detector regions excluding point sources identified using wavdetect as well as the galaxy cluster itself, and we used the CIAO tool deflare to identify periods of flaring. The resulting useful time on source, known as the “good time interval” (GTI), is reported in Table 4. We then extracted new events files using those GTIs, and those clean event files were used for all further X-ray analysis.

For the X-ray spectral analyses used to produce $T_{\textrm{\footnotesize e}}$ maps, the stowed ACIS background files were rescaled to match the high energy (10–12 keV) count rates off source (again, excluding regions with point source and galaxy cluster emission). These rescaled backgrounds were used as backgrounds in the spectral analysis. The regions used for spectroscopy were selected using the contour binning method of Sanders (2006).666Known as contbin, https://github.com/jeremysanders/contbin. The parameters were chosen to ensure each region had sufficient counts (typically $\gtrsim 3000$ background-subtracted counts from the inner portion of the galaxy cluster, though MACS J2129.4$-$0741 had $\sim 1800$) per spectral bin for reliable spectroscopy. The spectral analysis was carried out jointly for all available datasets in Sherpa (Freeman et al., 2001), using the xsmekal implementation of the Mewe-Kaastra-Liedahl (MeKaL) model. The hydrogen column density $N_{H}$ was fixed to the value found using the CIAO tool prop colden to obtain an interpolation of the Dickey & Lockman (1990) value at the galaxy cluster location. The redshift and abundance were also fixed in the analysis. In the case of abundance, several fits with abundance left free were also tested, and found not to differ significantly from fixing it to $Z=0.3~{}Z_{\odot}$.

4.4 HST

The HST images used to construct the lens models were already reduced, typically using standard procedures (most notably multidrizzle, see Koekemoer et al. 2011). For all of the galaxy clusters, previous lensing analyses exist, including multiple image constraints. The galaxy clusters were modeled here using parametrized forms, namely, double pseudo isothermal elliptical mass distributions for the galaxy cluster galaxies following common scaling relations, and elliptical NFW haloes for the galaxy cluster dark matter clumps. For the CLASH galaxy clusters, we adopted the Zitrin et al. (2015) “PIEMDeNFW” mass models. For the HST Frontier Fields galaxy cluster MACS J0717.5$+$3745 we remade and updated the model that is available on the HST Frontier Fields website.777https://frontierfields.org For modeling the RELICS galaxy clusters, we adopted the constraints from Cibirka et al. (2018) and Acebron et al. (2019) and we constructed a model for MACS J0018.5$+$1626 based on the constraints identified by Zitrin et al. (2011). For the remaining two galaxy clusters, we constructed models based on the multiple-image constraints listed in Richard et al. (2010) and Zitrin et al. (2011). Then, as our aim here was to supply maps to lens the CIB at radii well beyond the strong-lensing regime, and since our models were constructed from analytic, parametric forms, we then regenerated the strong-lensing models using the best-fit parameters from the above, but covering a larger field of view extending to the weak lensing regime. It should therefore be noted that these models have been extrapolated, as they were only constrained using data from within the HST field of view (solely strong lensing constraints, except for the CLASH galaxy clusters where HST weak lensing constraints were also used, see Zitrin et al. 2015). We regenerated all of the models onto a $16\arcmin\times 16\arcmin$ map, adopting a resolution of 0.25″ per pixel. Using these extended lens models we ray-traced different realizations of the background DSFGs that comprise the CIB, as detailed in Section 4.1.

5 Bulk Galaxy Cluster Velocities

5.1 Method and Results

Using the images produced in Section 4.1, from which radio galaxies, DSFGs, and the average lensing-induced CIB signal deficit were subtracted, we fitted a parametric model of the SZ effect signal to the data. First, an elliptical gNFW model, with power law exponents $\alpha$, $\beta$, and $\gamma$ fixed to the values found by Arnaud et al. (2010), was simultaneously fitted to the 140 GHz, 270 GHz, and Planck $y$–map data assuming a purely thermal SZ effect spectrum (i.e., zero kinematic SZ effect signal). As in Section 4.1, the image transfer functions and PSFs were fully accounted for in this fit. After this initial fit, which was used to determine the two-dimensional shape of the SZ effect signal, we then performed additional fits, separately to the 140 GHz and 270 GHz data, where only the normalization of the gNFW model was allowed to vary. This normalization was expressed in terms of the average surface brightness, in MJy sr*-1*, within an aperture centered on the galaxy cluster and extending to a radius of $R_{2500}$.

After adding the best-fit SZ effect model to each of the 1000 noise realizations for each galaxy cluster, an analogous two-step fit was performed, and the spread of normalization values obtained from these 1000 fits was used to estimate the uncertainty on that parameter. In addition, the $\chi^{2}$ values obtained from these 1000 fits were used to empirically determine the fit quality based on a probability to exceed (PTE) using the procedure described in Sayers et al. (2011). The average PTE for the 10 clusters was 0.37, and only one cluster had a PTE below 0.19 (RX J1226.9$+$3332, which had a PTE of 0.01). Therefore, even though many of these clusters are complicated mergers, the elliptical gNFW model was sufficient to describe our data given their noise and angular resolution. Furthermore, unlike X-ray observations, which are proportional to ICM density squared, the SZ effect data are linearly proportional to the ICM parameters (pressure for the thermal SZ effect and LOS velocity weighted by electron number density for the kinematic SZ effect). As a result, merger-induced ICM sub-structures, which can significantly bias similar bulk fits of smooth models to X-ray data, are much less problematic for fits to SZ effect data (e.g., Motl et al. 2005; Kay et al. 2012).

The SZ effect brightness values obtained from the above procedure were then used to constrain the overall bulk velocity of each galaxy cluster via the kinematic SZ effect. Specifically, we assumed that each galaxy cluster was moving with a single bulk LOS velocity and that its ICM was isothermal, with an electron temperature $T_{\textrm{\footnotesize e}}$ equal to the spectroscopic X-ray temperature measured by Chandra within $R_{2500}$. Given the assumption of an isothermal ICM with $T_{\textrm{\footnotesize e}}$ measured by Chandra, the total brightness from the thermal and kinematic SZ effect signals could be completely specified in terms of the electron optical depth $\tau_{\textrm{\footnotesize e}}$ and the bulk LOS velocity $v_{z}$. We used the SZpack software described in Chluba et al. (2012, 2013) to compute the SZ effect brightness for a given set of parameters, including relativistic corrections and assuming the effective thermal and kinematic SZ effect band centers given in Table 2. The results of these fits are shown in Figure 3 and summarized in Table 7.

The typical per-cluster uncertainty on the value of $v_{z}$ we obtained from these fits was 500–1000 km s*-1*, which was a factor of 2–4 larger than the typical expected $v_{z}$ (e.g., Evrard et al. 2002; Hernández-Monteagudo & Sunyaev 2010; Nagai et al. 2013). Not surprisingly, given these uncertainties, we did not detect a significant non-zero value of $v_{z}$ for any single galaxy cluster. To characterize the galaxy cluster ensemble as a whole, we computed the inverse variance weighted sample mean \langle$$v_{z}$$\rangle=430\pm 210 km s*-1*. However, this simple calculation did not account for the intrinsic cosmological variation in the value of $v_{z}$, and so we also computed the sample average velocity using a more sophisticated fit based on the linmix_err formalism of Kelly (2007). From these fits, we obtained a sample average velocity of \langle$$v_{z}$$\rangle=460\pm 300 km s*-1* and an intrinsic scatter with an rms of $\sigma_{\textrm{\small int}}=470\pm 340$ km s*-1*.

As expected, the mean velocity we obtained for our sample from both methods was consistent with zero, although the weighted mean differed at a significance of $\simeq 2\sigma$. While our uncertainty on the mean velocity was better than the pioneering measurements from SuZIE (Benson et al., 2003) and the value of $\pm 383$ km s*-1* obtained by Lindner et al. (2015) for a similar analysis of eleven galaxy clusters using data from the ACT and LABOCA, it was notably larger than the value of $\pm 60$ km s*-1* obtained from a Planck analysis of $\sim 1750$ X-ray-selected galaxy clusters (Planck Collaboration et al., 2014).

Our best-fit value for the intrinsic cluster-to-cluster scatter was consistent with the simulation-based expectation of $\sim 250$ km s*-1* (e.g., Evrard et al. 2002; Hernández-Monteagudo & Sunyaev 2010; Nagai et al. 2013), although with a somewhat large uncertainty of $\pm 340$ km s*-1*. However, we note that this uncertainty was comparable to what was obtained from Planck-based analyses of large samples of X-ray-selected galaxy clusters (i.e., $<800$ km s*-1* at 95% confidence in Planck Collaboration et al. 2014 and $350\pm 270$ km s*-1* in Planck Collaboration et al. 2018) and slightly better than the upper limit of 1450 km s*-1* obtained by Lindner et al. (2015).

5.2 Potential Sources of Bias

We note that the value of $v_{z}$ we obtained represents the average LOS velocity within $R_{2500}$. However, internal velocities in the ICM are expected to be comparable to the overall galaxy cluster peculiar velocity, even when the galaxy cluster is relatively relaxed. On average, these internal motions were not expected to produce a bias in the measured value of $v_{z}$, although they were expected to introduce an rms dispersion of $\simeq 50$–100 km s*-1*, depending on the orientation and the dynamical state of the galaxy cluster (Nagai et al., 2003). This dispersion is roughly one order of magnitude below our typical measurement uncertainty per cluster, and was therefore not included in our analysis.

The galaxy clusters in our sample are not isothermal, and so, in general, our assumption of an isothermal ICM produced some slight biases in our results (see, e.g., Chluba et al. 2013). Because the relativistic corrections to the SZ effect signal are non-linear with respect to $T_{\textrm{\footnotesize e}}$, the signal from an isothermal galaxy cluster will not in general be equal to the signal from a non-isothermal galaxy cluster with the same mean $T_{\textrm{\footnotesize e}}$. To estimate the potential bias from this effect, we computed the the expected SZ effect signal within $R_{2500}$ for the least isothermal galaxy cluster in our sample, MACS J0717.5$+$3745, using both the isothermal assumption and the 34 different values of $T_{\textrm{\footnotesize e}}$ within the separate contbin regions for that cluster. Even with $T_{\textrm{\footnotesize e}}$ ranging from 2 to 24 keV within those separate contbin regions, the fractional difference between the SZ effect signals computed using the two methods was only 0.2% at 140 GHz and 0.7% at 270 GHz. Assuming a similar $T_{\textrm{\footnotesize e}}$ structure along the LOS, this calculation indicates that the potential bias from our isothermal assumption was $\lesssim 1$% for all of the galaxy clusters in our sample.

Another, potentially larger source of bias was due to our use of X-ray spectroscopy from Chandra to determine the values of $T_{\textrm{\footnotesize e}}$. We note that the thermal SZ effect signal, and relativistic corrections to the SZ effect signals, depend on the LOS mass-weighted value of $T_{\textrm{\footnotesize e}}$. Within $R_{2500}$, hydrodynamical simulations indicate that the value of $T_{\textrm{\footnotesize e}}$ inferred from fitting an X-ray spectrum with a thermal emission model typically differ from the LOS mass-weighted $T_{\textrm{\footnotesize e}}$ at the level of 4–7% (Nagai et al. 2007b; see also Rasia et al. 2014). We did not attempt to correct for this difference in our analysis, although we note that it was sub-dominant compared to our assumed X-ray calibration uncertainty of 10%.

In addition, the well established difference in calibration between the two great X-ray observatories, Chandra and XMM-Newton, may also suggest a potential bias in our results. For the clusters in our sample, with $T_{\textrm{\footnotesize e}}$ generally between 5–10 keV, Chandra has been shown to systematically measure $T_{\textrm{\footnotesize e}}$ values $\simeq 10$–20% higher than XMM-Newton (e.g., Reese et al. 2010; Mahdavi et al. 2013; Donahue et al. 2014; Schellenberger et al. 2015; Madsen et al. 2017). While it is not clear which observatory has the more accurate calibration, this difference implies calibration uncertainties that may exceed the 10% rms we assumed in our analysis. Reconciling the Chandra/XMM-Newton calibration was beyond the scope of this work, but a relatively accurate post facto correction can be applied to our results if future work is able to better determine the effective area of Chandra. Because the relativistic corrections to the SZ effect signals were relatively small for our data (e.g., $\sim 10$% changes in $T_{\textrm{\footnotesize e}}$ would result in $\sim 1$% changes to the relativistic corrections), the spectral shapes of the thermal and kinematic SZ effect signals will remain nearly identical for small changes in $T_{\textrm{\footnotesize e}}$. Therefore, the thermal and kinematic SZ effect brightnesses obtained from our analysis would remain largely unchanged. As a result, if the value of $T_{\textrm{\footnotesize e}}$ changes by a factor of $1+\delta_{T}$, then, to good approximation, the value of $v_{z}$ will also change by a factor $1+\delta_{T}$ and the value $\tau_{\textrm{\footnotesize e}}$ will change by a factor of $1/(1+\delta_{T})$.

6 Resolved SZ Effect Imaging

As detailed in Czakon et al. (2015), it is possible to deconvolve the filtering effects described in Section 4.1 to obtain an unbiased image of the galaxy cluster SZ effect signal. One subtlety is that the filtering completely removes the mean signal level of the image, and so it must be determined using an independent measurement. In general agreement with the procedure of Czakon et al. (2015), we used the elliptical gNFW fits from Section 5 in order to determine the mean signal level of the unfiltered images. However, one important difference in this work was the addition of Planck $y$–map data in constraining the gNFW fits, as it was far more sensitive to the large angular scale SZ effect signal than the Bolocam data. Specifically, for this analysis we added a constant signal separately to the 140 and 270 GHz unfiltered images such that the average surface brightness within $R_{2500}$ was equal to the value obtained from the gNFW model fit. Example images are shown in Figure 5.

After we obtained these mean-corrected unbiased images, we then convolved them with a Gaussian kernel to obtain a common resolution of 70″ FWHM. While it would have been possible to use a resolution of 59″ FWHM, 70″ was chosen as a reasonable compromise between retaining spatial fidelity and filtering noise on small angular scales. From these images, we then fitted an SZ effect spectrum to each map pixel using the same procedure applied to the bulk galaxy cluster fits that were described in Section 5. Resolved maps of $T_{\textrm{\footnotesize e}}$ using Chandra X-ray spectroscopy were used to estimate $T_{\textrm{\footnotesize e}}$ within each pixel. From these fits, we then reconstructed resolved images of the thermal and kinematic SZ effect signals, which were then combined with the $T_{\textrm{\footnotesize e}}$ map to obtain images of the electron optical depth $\tau_{\textrm{\footnotesize e}}$ (see Figure 6) and the LOS velocity $v_{z}$ (see Figure 7).

For all ten galaxy clusters in our sample, the optical depth was imaged at high significance, with a peak S/N of more than 5. However, the most significant excursion identified in any of the velocity images had a S/N of 3, and it was coincident with the merging sub-cluster in MACS J0717.5$+$3745 previous described in Sayers et al. (2013a). Therefore, in nine of the ten clusters we were unable to detect the LOS velocity of any single sub-structure. To further search for evidence of underlying LOS velocity sub-structure below our detection limit within any single resolution element, we also computed the rms of the $v_{z}$ map over the galaxy cluster face within $R_{2500}$, $\sigma_{map}$. We then computed an identical rms from each of the 1000 noise realizations for each cluster ($\sigma_{noise}$), which provided an estimate of the expected rms in the absence of any underlying LOS velocity variations. We estimated the true internal $v_{z}$ rms as the difference between the measured rms and the expected rms due to noise (i.e., $\sigma^{2}_{vz}=\sigma^{2}_{map}-\langle\sigma_{noise}\rangle^{2}$). The distribution of $\sigma_{noise}$ values was also used to empirically determine confidence regions for the value of $\sigma_{vz}$. The resulting constraints on the rms of $v_{z}$ within $R_{2500}$ for each galaxy cluster are given in Table 7 and plotted in Figure 8.

All four of the clusters previously identified as likely LOS mergers had a non-zero measured $\sigma_{vz}$ (at a significance of $\simeq 2\sigma$ for Abell 0697, MACS J0018.5$+$1626, and MACS J2129.4$-$0741 and at a significance of $\simeq 4\sigma$ for MACS J0717.5$+$3745). The inferred $v_{z}$ rms for these clusters was $\gtrsim 1000$ km s*-1*, $\simeq 3$ times higher than expected from simulations of similar mass clusters (e.g., Nagai et al. 2013). While of modest statistical significance, our measurements were therefore consistent with a scenario where each of these four clusters is undergoing a merger along the LOS, which would boost the value of $\sigma_{vz}$. In contrast, the six clusters previous identified as likely POS mergers or relaxed all had a measured $\sigma_{vz}$ consistent with zero. At a confidence level of 95%, the $v_{z}$ rms for these clusters was $\lesssim 1000$–1500 km s*-1*. Based on the previously inferred merger geometry for the ten clusters in our sample, our SZ effect measurements were therefore able to distinguish LOS mergers from POS mergers and relaxed clusters.

One of the galaxy clusters in our sample, MACS J0717.5$+$3745, has been the target of several previous kinematic SZ effect studies, most notably by Sayers et al. (2013a) and Adam et al. (2017). Sayers et al. (2013a) used nearly identical data to those used in our study, although they included X-ray observations from XMM-Newton and they did not use Planck SZ effect data. They also used a much more individualized SZ effect analysis based on a spatial template derived from the X-ray data and a focus solely on the signal within 60″ diameter apertures centered on sub-clusters “B” and “C”. Their “direct integration” results are the most comparable to those obtained in our more general SZ effect analysis, and they obtained best-fit $v_{z}$ values of $+2550\pm 1050$ km s*-1* towards “B” and $-500\pm 1600$ km s*-1* towards “C”. At the same positions in our $v_{z}$ map, we obtained values of $2100\pm 700$ km s*-1* and $-400\pm 800$ km s*-1*. The shift to a smaller positive $v_{z}$ for sub-cluster “B” was driven mainly by our correction for the lensing-induced deficit in the CIB, which was not included in the analysis of Sayers et al. (2013a). This also drove a shift towards a larger negative $v_{z}$ for “C”, but this shift was more than compensated for by the significantly lower $T_{e}$ obtained in our analysis, which resulted in a larger best-fit $\tau_{\textrm{\footnotesize e}}$ and subsequently smaller magnitude for $v_{z}$. The smaller uncertainties obtained in our analysis were due to the combination of: a larger aperture (70″ compared to 60″); larger best-fit values for $\tau_{\textrm{\footnotesize e}}$, particularly for “C”; improved calibration; and, most significantly, the inclusion of Planck SZ effect data to better constrain the large angular scale signal.

Using completely independent SZ effect measurements, Adam et al. (2017) measured best-fit $v_{z}$ values of $+6600^{+3200}_{-2400}$ km s*-1* and $-4100^{+1600}_{-1100}$ km s*-1* for sub-clusters “B” and “C”. These values are in modest ($\lesssim 2\sigma$) tension with the values we obtained in our analysis, although we note that Adam et al. (2017) found $v_{z}$ equal to $+2100^{+500}_{-450}$ km s*-1* for sub-cluster “B” using an alternate analysis which included stronger X-ray priors, fully consistent with our measurement. Furthermore, the NIKA SZ effect observations used by Adam et al. (2017) had a factor of $\simeq 3$ finer angular resolution compared to our Bolocam/AzTEC data, better isolating the sub-clusters and producing more significant excursions in their $v_{z}$ map of the galaxy cluster ($5.1\sigma$ and $3.4\sigma$ for “B” and “C”).

7 Summary

We have used observations from Bolocam and AzTEC to image the SZ effect signal towards a sample of ten galaxy clusters at 140 and 270 GHz. In support of these data, we have also made use of a number of additional observations. The Planck all-sky $y$–maps were used to help constrain the large-angular scale signal in order to obtain spatial templates of the SZ effect signal. In addition, three-band Herschel–SPIRE imaging was used to subtract the emission from DSFGs, which was significant compared to the SZ effect signal at 270 GHz. Furthermore, HST data were used to obtain detailed mass models for each galaxy cluster in order to properly account for lensing of the background CIB. Finally, Chandra X-ray spectroscopic imaging was used to obtain resolved maps of the ICM temperature $T_{\textrm{\footnotesize e}}$.

From this analysis, we produced galaxy cluster-averaged fits to the SZ effect brightness at 140 and 270 GHz in order to constrain the average optical depth and bulk LOS velocity $v_{z}$ within $R_{2500}$. Our typical measurement uncertainties on $v_{z}$ were 500–1000 km s*-1*, a factor of 2–4 larger than the typical values of $v_{z}$ expected from simulations. We did not detect $v_{z}$ at high significance in any single galaxy cluster, and the ensemble average velocity was consistent with zero, particularly when intrinsic scatter was accounted for in the fit. When fitting for the intrinsic scatter, we did not obtain a significant detection, but we did find an upper limit competitive with those produced by statistical stacks in CMB survey data.

In addition to fitting for the galaxy cluster-average SZ effect brightness, we also produced images of the electron optical depth $\tau_{\textrm{\footnotesize e}}$ and the LOS velocity $v_{z}$ with a resolution of 70″. In all cases, $\tau_{\textrm{\footnotesize e}}$ was detected at high significance near the galaxy cluster center. We did not obtain a significant detection of $v_{z}$ within any single resolution element for any of the galaxy clusters in our sample, with the exception of the previously identified sub-component of MACS J0717.5$+$3745. However, all four of the clusters previously identified as likely LOS mergers showed a $v_{z}$ rms greater than zero at a significance of $\gtrsim 2\sigma$, with $\sigma_{vz}\gtrsim 1000$ km s*-1* for these objects. This is a factor of $\simeq 3$ above the $v_{z}$ rms expected from simulated clusters of similar masses (e.g., Nagai et al. 2013), strongly indicating a boosted $\sigma_{vz}$ due to a LOS merger. In contrast, all six of the clusters previously identified as likely POS mergers or relaxed had $\sigma_{vz}$ consistent with zero and $\sigma_{vz}\lesssim 1000$–1500 km s*-1* at a 95% confidence level. Based on the previous characterizations of the merger geometries for these galaxy clusters, our SZ effect data were therefore able to distinguish between LOS mergers from POS mergers and relaxed clusters.

In addition to the ICM constraints obtained in our analysis, we also quantified the potential bias in measuring the SZ effect signal due to lensing of the background DSFGs that comprise the CIB. When individual bright DSFGs are identified and subtracted, lensing produces an on-average deficit in the surface brightness of the CIB. For the galaxy clusters in our sample, the total surface brightness of this deficit was typically $\simeq 15$% of the total surface brightness of the SZ effect signal at 270 GHz, although it was as large as 25% for one galaxy cluster (MACS J2129.4$-$0741).

In contrast to some other recent kinematic SZ effect analyses (e.g., Sayers et al. 2013a and Adam et al. 2017), we did not make use of X-ray data to model the shape of the ICM pressure or density, although we did make use of resolved temperature maps from spectroscopic Chandra X-ray observations. While such X-ray density and pressure information can be useful in breaking degeneracies in kinematic SZ effect measurements (e.g., Flender et al. 2017), they can be difficult to include for detailed studies of individual galaxy clusters with complicated merger geometries, as was the case for nine of the ten objects in our sample. For example, deprojections, which assume a spherical geometry, can only be applied in special cases for merging clusters (e.g., it was only used for one ICM sub-component in the analysis of Adam et al. 2017). Sayers et al. (2013a) avoided this complication by using the SZ effect data to constrain the LOS extent of an X-ray derived pseudo-pressure map of MACS J0717.5$+$3745, although the resulting constraints were only marginally better than those obtained in our current analysis without an X-ray template. We therefore decided for this work not to pursue an analysis based on X-ray maps of the ICM density or pressure.

Looking to the future, instruments like TolTEC (Bryan et al., 2018) will be able to provide much deeper SZ effect observations, at finer angular resolution, and in more observing bands. Scheduled to be installed in 2019 on the 50 meter Large Millimeter Telescope Alfonso Serrano (LMT) in México, TolTEC will simultaneously observe at 150, 220, and 280 GHz, providing images of the thermal and kinematic SZ effect signals while robustly detecting (and subtracting) the contaminating signals from DSFGs. Compared to the data used for this work, TolTEC promises an order of magnitude or more improvement in achievable map depth, opening up the prospect of high significance imaging of ICM velocity structures at $\lesssim 10$″ resolution. Longer-term concepts, such as AtLAST (Bertoldi, 2018) and CSST (Golwala, 2018), promise to deliver much larger fields of view than TolTEC and expanded spectral coverage (e.g., 90–400 GHz), further enhancing the potential scientific reach of detailed SZ effect studies.

8 Acknowledgments

Sunil Golwala, Erik Reese, Jack Sayers, and Grant Wilson were supported by NASA/NNX15AV66G. This work was partially supported by Consejo Nacional de Ciencia y Technología (CONACYT) project FDC-2016-1848. The scientific results reported in this article are based in part on data obtained from the Chandra Data Archive. This research has made use of software provided by the Chandra X-ray Center (CXC) in the application packages CIAO, ChiPS, and Sherpa. We thank Georgiana Ogrean for providing her X-ray reprocessing and analysis scripts. We thank the referee for providing several useful suggestions that have improved this manuscript.