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

All of the galaxy clusters in our sample were imaged with Bolocam at 140 GHz,¹³ and all of those data have been used in previous analyses (e.g., Sayers et al. 2013c; Czakon et al. 2015) and are publicly available.¹⁴ 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 592. The data were collected in 10 minute observations using a sinusoidal Lissajous pattern with differing periods in the R.A. and decl. 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 ≃5'', 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 ∼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 ≃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).

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 332. 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%.

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 had a solid angle that corresponded to a Gaussian with a FWHM of 304. 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).

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 (Egami et al. 2010). Herschel–SPIRE was a three-band photometric imager operating at 600, 850, and 1200 GHz with PSFs having FWHMs of 181, 252, and 366 (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.

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' × 8', observations with either ACIS-I or ACIS-S covered a sufficiently large field of view for this analysis.

We reconstructed lens models for the 10 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; A0697, 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 A1835 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 005–006 per pixel. Details of the lens modeling are given in Section 4.4.

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 ≃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. 2013c).

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. (2016b). 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 Navarro–Frenk–White (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. (2016a) 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 α, β, and γ to the best-fit values of Arnaud et al. (2010). For the radial scales typically probed by our data, ≃0.3R₂₅₀₀–3.0R₂₅₀₀,¹⁵ this model had sufficient freedom to provide a good fit quality (see Sayers et al. 2011; 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 α, β, and γ (e.g., Planck Collaboration et al. 2013; Sayers et al. 2013c; 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 signal-to-noise ratio (S/N) > 4 from the resulting images. We typically detected ≃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. (2013b) 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. (2013b) 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 graybody spectral energy distribution (SED) of the form

$$\begin{eqnarray}&&F(\nu )={F}_{0}(1-{e}^{(-\nu /{\nu }_{0}{)}^{\beta }}\,)B(\nu ,{T}_{{\rm{d}}}),\end{eqnarray} \tag{ 1 }$$

where the values of the normalization F₀ and the dust temperature ${T}_{{\rm{d}}}$ were varied, ν is the observed frequency, ν₀ = 3000 GHz (Draine 2006), β = 1.95 is the dust emissivity spectral index,¹⁶ and $B(\nu ,{T}_{{\rm{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 graybody 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 graybody 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. (2013c). Since the quality of the data used in this work was nearly identical to that used by Sayers et al. (2013c), the same overall implications were also true and are summarized here. In particular, all sources brighter than ≃1 mJy at 270 GHz were detected, and some sources were detected down to a limit of ≃0.1 mJy at 270 GHz. In aggregate, these detected sources represent ≃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 ≃5%–10% of the total CIB. Even after subtracting ≃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 ≃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 intrinsically 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 ≃15% of the SZ effect brightness at 270 GHz (see Table 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 (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 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 an S/N threshold on those data. This typically resulted in ≃100 DSFG candidates per galaxy cluster (see Table 5).

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.¹⁷ 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}_{{\rm{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).¹⁸ 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 ∼1800) per spectral bin for reliable spectroscopy. The spectral analysis was carried out jointly for all available data sets in Sherpa (Freeman et al. 2001), using the xsmekal implementation of the Mewe–Kaastra–Liedahl 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_⊙.

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 parameterized 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, which is available on the HST Frontier Fields website.¹⁹ 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' × 16' map, adopting a resolution of 025 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.

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 α, β, and γ 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⁻¹, within an aperture centered on the galaxy cluster and extending to a radius of R₂₅₀₀ (see Figure 2).

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 χ² 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}_{{\rm{e}}}$ equal to the spectroscopic X-ray temperature measured by Chandra within R₂₅₀₀. Given the assumption of an isothermal ICM with ${T}_{{\rm{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 }_{{\rm{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 Figures 3 and 4 and summarized in Table 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 7.  Derived ICM Parameters

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

Note. Best-fit values and 68% confidence intervals for the ICM parameters derived in our analysis. The temperature constraints were obtained from Chandra (not including the assumed 10% systematic uncertainty), and the optical depth and bulk velocity constraints were obtained from our SZ effect fits. The internal v_z rms constraints were obtained from the resolved SZ effect maps within R₂₅₀₀. The 68% confidence interval for six clusters is consistent with an internal v_z rms of zero, and 95% confidence level upper limits are given for these clusters.

Download table as:  ASCIITypeset image

The typical per-cluster uncertainty on the value of v_z we obtained from these fits was 500–1000 km s⁻¹, 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⁻¹. 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⁻¹ and an intrinsic scatter with an rms of σ_int = 470 ± 340 km s⁻¹.

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 ≃2σ. While our uncertainty on the mean velocity was better than the pioneering measurements from SuZIE (Benson et al. 2003) and the value of ±383 km s⁻¹ obtained by Lindner et al. (2015) for a similar analysis of 11 galaxy clusters using data from the ACT and LABOCA, it was notably larger than the value of ±60 km s⁻¹ obtained from a Planck analysis of ∼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 ∼250 km s⁻¹ (e.g., Evrard et al. 2002; Hernández-Monteagudo & Sunyaev 2010; Nagai et al. 2013), although with a somewhat large uncertainty of ±340 km s⁻¹. 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., $\lt 800$ km s⁻¹ at 95% confidence in Planck Collaboration et al. 2014 and 350 ± 270 km s⁻¹ in Planck Collaboration et al. 2018) and slightly better than the upper limit of 1450 km s⁻¹ obtained by Lindner et al. (2015).

We note that the value of v_z we obtained represents the average LOS velocity within R₂₅₀₀. 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 ≃50–100 km s⁻¹, 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 nonlinear with respect to ${T}_{{\rm{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}_{{\rm{e}}}$. To estimate the potential bias from this effect, we computed the the expected SZ effect signal within R₂₅₀₀ 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}_{{\rm{e}}}$ within the separate contbin regions for that cluster. Even with ${T}_{{\rm{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}_{{\rm{e}}}$ structure along the LOS, this calculation indicates that the potential bias from our isothermal assumption was ≲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}_{{\rm{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}_{{\rm{e}}}$. Within R₂₅₀₀, hydrodynamical simulations indicate that the value of ${T}_{{\rm{e}}}$ inferred from fitting an X-ray spectrum with a thermal emission model typically differs from the LOS mass-weighted ${T}_{{\rm{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}_{{\rm{e}}}$ generally between 5 and 10 keV, Chandra has been shown to systematically measure ${T}_{{\rm{e}}}$ values ≃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., ∼10% changes in ${T}_{{\rm{e}}}$ would result in ∼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}_{{\rm{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}_{{\rm{e}}}$ changes by a factor of 1 + δ_T then, to a good approximation, the value of v_z will also change by a factor 1 + δ_T and the value ${\tau }_{{\rm{e}}}$ will change by a factor of 1/(1 + δ_T).