MUSTANG HIGH ANGULAR RESOLUTION SUNYAEV–ZEL'DOVICH EFFECT IMAGING OF SUBSTRUCTURE IN FOUR GALAXY CLUSTERS

P. M. Korngut; S. R. Dicker; E. D. Reese; B. S. Mason; M. J. Devlin; T. Mroczkowski; C. L. Sarazin; M. Sun; J. Sievers

doi:10.1088/0004-637X/734/1/10

1. INTRODUCTION

The Sunyaev–Zel'dovich effect (SZE) in clusters of galaxies arises from inverse Compton scattering of cosmic microwave background (CMB) photons off hot electrons in the intracluster medium (ICM; Sunyaev & Zel'dovich 1972). The magnitude of this effect is redshift independent and directly proportional to the line-of-sight integrated pressure of the plasma. At frequencies ≲ 218 GHz, it is manifested as a decrement in CMB intensity. Over the past two decades, measurements of the SZE in clusters of galaxies have been used to probe a wide range of cosmological and astrophysical questions. It has been used by dedicated surveys to search for clusters (e.g., Hincks et al. 2010; Vanderlinde et al. 2010; Menanteau et al. 2010; Marriage et al. 2010), combined with X-ray data to measure the Hubble flow (e.g., Mason et al. 2001; Reese et al. 2002; Bonamente et al. 2006) and to derive physical cluster properties from radial profiles (e.g., LaRoque et al. 2003; Bonamente et al. 2006; Mroczkowski et al. 2009). For reviews of the SZE and its applications, see Birkinshaw (1999) and Carlstrom et al. (2002).

Measurements of the SZE at high angular resolution are difficult because of the large collecting areas required. Nearly all measurements currently in the literature have effective angular resolution ≳ 1'. These large angular scales (corresponding to ∼365 kpc at z = 0.5) are useful for measuring cluster properties within large clustercentric radii (e.g., Nord et al. 2009) but are unable to resolve cluster substructure within the ICM.

High-resolution X-ray imaging from Chandra and XMM-Newton in the last decade opened a new window to cluster physics. Objects once thought to be spherically symmetric and relaxed have been shown to display evidence of interesting phenomena which provide insight to the complicated dynamics at play in these structures. Among these are shocks and cold fronts induced by recent mergers (e.g., Markevitch & Vikhlinin 2007), cavities and heating caused by active galactic nucleus (AGN) interactions (e.g., McNamara et al. 2005), and sharp surface brightness edges caused by gas sloshing (e.g., ZuHone et al. 2010).

High-resolution images of the SZE in clusters provide a new tool which, when combined with X-ray measurements, can constrain complicated physics in galaxy clusters. This is particularly true of the high-redshift universe as the X-ray surface brightness data (proportional to the product of the density squared and square root of temperature integrated along the line of sight) suffer from cosmological dimming.

The potential of resolved SZE was first demonstrated by Komatsu et al. (2001), who used the Nobeyama 45 m to image RXJ1347−1145, a massive X-ray luminous cluster previously thought to be relaxed and spherically symmetric (Schindler et al. 1997). The asymmetry revealed by their work was the first indication that the system was disturbed, and it is now believed to have undergone a recent merger (Kitayama et al. 2004). This has since been confirmed in the SZE by MUSTANG (Mason et al. 2010). More recently, other groups have made high-resolution SZE images in CL J0152−1347 (Massardi et al. 2010, at ∼35'' resolution) and the Bullet Cluster (Savyasachi Malu et al. 2010, at ∼30'' resolution).

In this work, we present measurements taken with the Multiplexed SQUID/TES Array at Ninety Gigahertz (MUSTANG) receiver on the 100 m Robert C. Byrd Green Bank Telescope (GBT). The large collecting area of the GBT combined with the focal plane array of bolometers make this system ideal for probing substructure in clusters through the SZE.

The clusters MACS0744, RXJ1347, and CL1226 were selected because they displayed signs of merger activity in previous measurements and therefore were likely to contain small-scale features created by merger processes. In contrast, A1835 was specifically targeted because it was expected to be representative of a typical relaxed cluster. All uncertainties quoted in this paper are 68% confidence and we assume a cosmology where H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.73, and Ω_M = 0.27.

2. INSTRUMENT AND OBSERVATIONS

2.1. MUSTANG

MUSTANG is a focal plane camera with an 8 × 8 array of Transition Edge Sensor (TES) bolometers built for the Gregorian focus of the 100 m GBT. It has 18.4 GHz of bandwidth centered on 90 GHz. The array has a 0.63fλ pixel spacing which yields a well-sampled instantaneous field of view (FOV) of 42'' on the sky. More detailed information on the instrument can be found on the MUSTANG Web site⁶ and in Dicker et al. (2008, 2009).

2.2. Observations

Data presented here were obtained during the winter/spring of 2009 and 2010. The cluster signal was modulated predominantly in a "Lissajous daisy" scan pattern. This strategy was designed to move the telescope with high speed (∼0 farcm 5 s⁻¹) without drastic accelerations which can induce feed arm instabilities and pointing wobble. Faster scan rates move the sky signal to frequencies above the low-frequency (1/f) noise from the atmosphere and internal fluctuations. The GBT bore-sight trajectory during one of these scans is displayed in Figure 1. This observation pattern is centrally weighted and produces maps with radially increasing noise levels. While the particular scan shown in Figure 1 contains information in a ∼6' diameter region, only the central ∼2' is well covered. To improve the sky sampling in cluster cores, each object was mapped with five tiled pointing centers: one centered on the X-ray surface brightness peak and others offset to the north, south, east, and west by approximately one instantaneous FOV (∼42''). Other scan patterns with more uniform sky coverage such as the "billiard ball" scan described in Dicker et al. (2009), Cotton et al. (2009), and Mason et al. (2010) were used as well. This alternative strategy has the advantage of uniform noise across the map, but at the cost of slower telescope velocity with sharper turnarounds.

**Figure 1.** GBT bore-sight trajectory during a Lissajous daisy scan pattern. This observing strategy was used for the majority of the observations presented. The movement of the source across the sky during the scan has been subtracted. The shaded box indicates the instantaneous FOV of MUSTANG.
Download figure:
Standard image High-resolution image

At the start of each session, out-of-focus (OOF) holography was carried out using a bright (∼1 Jy) unresolved source. This technique, described in detail in Nikolic et al. (2007), consists of mapping a compact source with the GBT secondary in three positions relative to the primary: nominally in focus and 3λ on either side. An automated real-time analysis uses the measured beam patterns to fit for phase errors in the telescope aperture. Corrections to the active surface of the GBT primary are calculated and applied along with pointing and focus offsets.

After every two Lissajous daisy scans (∼30 minutes on source), the beam profile was measured using a nearby bright compact quasar. If significant ellipticity or gain decrease is detected in the periodic beam measurements, the OOF procedure is repeated. These beam maps are used in image reconstruction to track fluctuations in the telescope gain, atmosphere, and pointing offsets. The ∼30 minute calibration timescale was chosen as it is characteristic of the thermal time constant of the telescope. The sources used as secondary calibrators for each cluster along with the total on-source integration times are presented in Table 1.

Table 1. Observation Summary

Cluster	z_r	Time	Secondary Calibrator
		(hr)
A1835	0.25	3.5	1415+1320
RX J1347.5−1145	0.45	3.3	1337−1257
MACS J0744.8+3927	0.69	5.8	0824+3916
CL J1226.9+3332	0.89	4.5	1159+2914

Download table as: ASCII Typeset image

Planets were used for absolute flux calibration and were mapped at least once per night. The fluxes of these primary calibrators are taken from Weiland et al. (2011). Several times a night, off-source scans with the telescope at rest and the internal calibration lamp (CAL) firing with a 0.5 Hz square wave were taken. These are used in analysis to fit for the gain of each pixel. The absolute flux of the data is calibrated to an accuracy of 15%.

3. MUSTANG DATA REDUCTION

A custom imaging algorithm implemented in IDL is used to produce maps from the time-ordered bolometer data. The data are heavily filtered to remove atmospheric signal prior to map making. The process is outlined below.

1.
Gain inhomogeneities across the detector array are flat-fielded using the nearest CAL scan. These data are also used to identify and mask unresponsive pixels (typically 10–15 out of 64).
2.
A template of the atmospheric signal is estimated from low-frequency fluctuations that are highly correlated across the array. This is constructed from an average of the time streams from all accepted pixels. The model is then low-pass filtered in Fourier space to separate the astronomical signal on small spatial scales from the atmospheric template. This filtering requires a characteristic frequency based on the noise properties of the data. The template is then subtracted in the time domain from each pixel. The effectiveness of this filter relies on the assumption that the celestial signal is not common mode, which is valid only in the limit of compact sources. SZE signal from clusters can be smooth and often extend for many arcminutes on the sky and not well approximated in this assumption. It is therefore essential to simulate and quantify the angular transfer function of the imaging pipeline.
3.
A low-order polynomial is fit and subtracted from each time stream. This further removes the long timescale fluctuations in the data.
4.
The data all contain a coherent 1.411 Hz signal. This is produced by fluctuations in optical load on the detectors caused by the thermal cycle of the camera's main cryogenic refrigerator. It is well approximated by a sinusoid and is removed at this stage.
5.
A per-pixel high-pass filter is applied in Fourier space. This aggressive technique removes all low-frequency spatial modes from the data indiscriminately. A characteristic frequency is defined at this stage as well. This will further affect the angular scales present in the reconstructed image.
6.
Individual detector weights are computed based on the noise characteristics of each detector after the processing described above. This is used to create an effective exposure time for each pixel on the sky.
7.
The time stream data are then binned on a 2'' × 2'' grid in right ascension and declination.

The MUSTANG images presented in this work have been optimized for peak signal to noise on the compact features of the clusters. As described above, there are two selectable filter parameters used in the map maker, one for the common mode template and the other for the per-pixel high pass. These selected frequencies correspond to spatial scales on the sky through the speed at which the signal is modulated by the telescope scan (usually ∼0 farcm 5 s⁻¹). The optimal filter for each object depends on the intrinsic structure of the source as well as the noise properties of the scans used in each observation. To determine the optimal filter for each map, parameter space is explored systematically by mapping each object with varied degrees of filtering. A single optimal value for each parameter is assumed for the entire data set on each cluster.

The noise in each map is defined by the standard deviation of all pixels in an off-source region free from obvious signal. From extensive Monte Carlo simulation, we find that this calculation provides a good measurement of the noise in the map as a whole, provided that it is scaled by the square root of the difference in the map weights of the areas in question.

The necessary filtering steps described above result in an attenuation of flux in the recovered map. The magnitude of this attenuation depends strongly on angular scale. Typically, structure larger than the instantaneous FOV is attenuated, and by ∼100'' is less than half the amplitude of the signal within the instantaneous FOV (see Figure 2). The angular transfer function for each object mapped is calculated using the specific scans and filter parameters selected to produce the map. When estimating the flux in a map or quantitatively comparing observational data to model fits, it is essential to correct for the effects of this transfer function (as we discuss in Section 7). A detailed description of the calculation of the transfer function is described in Mason et al. (2010).

**Figure 2.** Effective transfer function, as a function of radial angular scale, after the filtering used to produce the MUSTANG map of CL1226 reported here (see Figure 10). Attenuation was calculated using simulated observations of 10 fake skies. Each sky was generated at random to include all angular scales larger than the MUSTANG beam, out to 2', and to have an rms of 1 mJy. The real telescope pointings for CL1226 were used to generate randomized time streams. Noise taken from off-source scans was added. These data were then analyzed using the same pipeline used to produce our map of CL1226. The one-dimensional transfer function plotted here is the ratio of azimuthally averaged power spectra of the reconstructed and input fake skies. The angles subtended by the beam and instantaneous FOV are labeled for clarity (vertical dash-dotted and dashed lines).
Download figure:
Standard image High-resolution image

4. CHANDRA DATA REDUCTION

Archival Chandra data are reduced using CIAO version 4.2 and calibration database 4.2.0. Starting with the level 1 events file, standard corrections are applied along with light curve filtering and other standard processing (for reduction details, see Reese et al. 2010). Images are made in full resolution (0 farcs 492 pixels) and exposure maps are computed at 1 keV. When merging data from separate observations, images and exposure maps from each data set are combined and a wavelet-based source detector is used on the combined image and exposure map to find and generate a list of potential point sources. The list is examined and adjusted by eye and used for our point-source mask. A summary of the archival Chandra data used in this paper is presented in Table 2.

Table 2. Archival Chandra Data

Cluster	Time	ObsIDs
	(ks)
A1835	222	495, 496, 6880,
		6881, 7370
MACS J0744.8+3927	90	3197, 3585, 6111
CL J1226.9+3332	74	3810, 5014, 932

Download table as: ASCII Typeset image

5. RESULTS

5.1. MACS J0744.8+3927 (z = 0.69)

This massive high-redshift system, found in the Massive Cluster Survey (MACS) of the all sky ROSAT data (Ebeling et al. 2001a), has appeared in several studies using X-ray and SZE data (e.g., LaRoque et al. 2003, 2006; Ebeling et al. 2007). Unlike the other clusters in our sample, targeted multi-wavelength studies of MACS0744 are scarce in the literature. Kartaltepe et al. (2008) include this object in a red sequence galaxy distribution study of a subsample of 12 MACS clusters. They note that understanding the assembly dynamics of this system is made difficult by its complex morphology, which includes some evidence of a dense core in the X-ray images, and an elongated, doubly peaked distribution of red sequence galaxies. Additionally, this cluster field contains no significant compact radio sources, as determined in a search of the Faint Images of the Radio Sky at Twenty-Centimeters (FIRST; White et al. 1997) and NRAO VLA Sky Survey (NVSS; Condon et al. 1998) catalogs at 1.4 GHz, as well as a literature search for SZA, OVRO, and BIMA sources detected in this cluster field at ≈30 GHz (LaRoque et al. 2006; Mroczkowski et al. 2009).

5.1.1. MUSTANG Data

The SZE map produced from 5.8 hr of MUSTANG data is shown in Figure 3. It consists of a kidney-shaped ridge ∼25'' long in the north–south direction. From east to west, the structure is roughly the width of our beam, and thus is not resolved in this direction. The curvature of this feature is well described empirically as an 80 deg sector of an ellipse with an axial ratio of 1.25, with the minor axis and center of the observed SZE being 12 deg south of west on the sky.

farcs — **Figure 3.** SZE and X-ray images of MACS0744. Left: MUSTANG+GBT SZE at 135 FWHM effective resolution after smoothing. Contours are multiples of 0.5σ starting at 3σ. Center: *Chandra* X-ray surface brightness in the cluster core. The image has been smoothed with a 15 Gaussian. Right: composite image of *Chandra* X-ray and MUSTANG SZE. Blue and green are identical data on different logarithmic color scales. Red shows the MUSTANG SZE data. The kidney-shaped ridge revealed by MUSTANG is aligned concentrically with a sharp surface brightness discontinuity in the *Chandra* map.
Download figure:
Standard image High-resolution image

5.1.2. Chandra Data

The Chandra image is shown beside the MUSTANG map in Figure 3. It was produced from nearly 90 ks of combined archival data merged from ObsIDs 3197, 3585, and 6111 and reduced with the method described in Section 4. The core of this cluster displays an asymmetric X-ray surface brightness morphology with a sharp discontinuity on the western edge. The concave side of the SZE peak identified by MUSTANG is aligned concentrically with the convex edge of the surface brightness discontinuity in the X-ray.

5.1.3. X-Ray Surface Brightness Shock Modeling

The combined SZE and X-ray image morphology presented in Figure 3 is suggestive of a system dominated by a merger-driven shock front. Arriving at this conclusion based on the existing relatively low signal-to-noise ratio (S/N) X-ray and SZE data alone would be quite tenuous; however, the kidney-shaped ridge seen by MUSTANG combined with the sharp edge seen by Chandra is difficult to explain without invoking a shock-heating mechanism. We proceed to model the system in the framework of a shock front through a complementary analysis of X-ray and SZE in the approach outlined below.

1.
The elliptical geometry and location of the shocked gas is approximated from the SZE data.
2.
This geometry is used to fit a two-dimensional X-ray surface brightness profile with a model consisting of three regions: a cold intact core bordered by a cold front, a shock-heated region bordered by a shock front, and a pre-shock region. These correspond to I, II, and III, respectively, in Figure 4.
3.
X-ray spectroscopy is performed in each region to obtain the plasma temperature.
4.
Three-dimensional density and pressure models are produced from the surface brightness and spectral fits.
5.
The pressure model is integrated along the line of sight to produce a two-dimensional Compton y_C map.
6.
A mock SZE image is constructed at the resolution and with the angular extent of the MUSTANG map, and the model and data are compared.

**Figure 4.** Geometry and regions used for elliptical profiles and X-ray spectroscopy on MACS0744 overlaid on the *Chandra* surface brightness image. Green contours are (− 4.5, −5.5)σ SZE decrement. The three regions correspond to the cool intact core (I), the shock-heated gas (II), and pre-shock region (III). One X-ray point source has been excised from the pre-shock region. The borders of the wedge indicate the azimuthal range used in producing radial profiles.
Download figure:
Standard image High-resolution image

We model the X-ray emissivity as a power law, ε∝r^−p, within each region assuming an ellipsoidal geometry with two axes in the plane of the sky and one along the line of sight (see Appendix A for details). The model has eight parameters in total, two characteristic radii, and a normalization and power-law index in each of three regions. We perform a Markov Chain Monte Carlo (MCMC) analysis using Poisson statistics for the X-ray data (for analysis and statistics details see, e.g., Reese et al. 2000, 2002; Bonamente et al. 2006). Each chain is run for a million iterations. Convergence and mixing are checked by running two chains and comparing them against one another (Gelman & Rubin 1992; Verde et al. 2003). The choice of burn-in period does not significantly affect the results but for concreteness we report results using a burn in of 10,000 iterations. The model fit is limited to a wedge subtending 80° and extending from 10'' to 40'' from the nominal center. This region corresponds to the region of interest suggested by the SZE and X-ray data as discussed in Section 5.1.2.

Initial attempts to model all eight parameters at once were unsuccessful due to low S/N in these small regions, with the chains showing poor convergence. To limit the number of free parameters, we implement chains to determine the discontinuity radii, R_s1 and R_s2, individually and then fix those radii. This entails using a single discontinuity model, which has five parameters, rather than eight. The inner discontinuity radius, R_s1, is determined with single discontinuity chains using the entire fitting region. The outer discontinuity radius, R_s2, is fit with a single discontinuity model limiting the fitting region to larger radii than R_s1.

With both discontinuity radii in hand, the double discontinuity model chains are run with fixed characteristic radii. This is enough of a reduction of parameter space to produce converged chains. Best fit and 68% confidence level uncertainties are shown in Table 3. In this table, the parameter f is defined to be the ratio of the normalization of a given region over the normalization in the cool intact core (region I). Because the radial dependence of the model follows a power law with an exponent less than zero, the amplitudes quoted here are normalized at the cold-front radius, R_s1 = 14 farcs 19, to avoid a singularity at the origin. Figure 5 shows the X-ray surface brightness profile within the fitting region along with the best-fit model.

**Figure 5.** Top left: *Chandra* X-ray surface brightness elliptical profile (points) and the best-fit analytical model in MACS0744 (red line). Bottom left: residual of data and model in the top left. Blue lines show the best-fit characteristic radii for the cold front and shock front, 1419 and 1923, respectively. Top right: intrinsic electron number density model produced from the surface brightness fit on the left. Bottom right: pressure model produced from the above density model and the temperatures derived from *Chandra* spectroscopy. The shaded regions show the uncertainty based on the spectroscopically measured temperatures. Radii in this figure are elliptical and follow the conventions described in Appendix A.
Download figure:
Standard image High-resolution image

Table 3. Best-fit Parameters for the Shock Model in MACS0744

Region	f	p	k_BT_e
			(keV)
I	1	0.913^+0.379_{− 0.285}	8.2^+1.6_{− 1.2}
II	0.480^+0.124_{− 0.084}	0.986^+0.559_{− 0.349}	19.7^+9.7_{− 5.9}
III	0.406^+0.086_{− 0.063}	1.151^+0.041_{− 0.040}	8.7^+1.1_{− 0.8}

Download table as: ASCII Typeset image

We also ran MCMC fits modeling a constant X-ray background in addition to the shock model. It has no statistically significant effect on the shock model results. This is not surprising as the X-ray background is over an order of magnitude down in surface brightness compared to the cluster signal at the outermost radius considered in the fit. The X-ray background becomes even less important toward the inner radii where the cluster signal rises.

To produce Compton y_C maps, the three-dimensional pressure model was numerically integrated along the line of sight using Equation (A17) out to an elliptical radius of 60'', where the single power-law model becomes a poor description of the X-ray data. This map is then used to produce a predicted SZE image at 90 GHz. After convolving with the GBT beam, the angular transfer function of the analysis pipeline is applied to the model in Fourier space and compared to the measured SZE data. Since the model is only valid in a specified range of angles about the center of the ellipse, the remaining sky was assumed to be well described by the double-β model of LaRoque et al. (2006) and a single temperature of 8.0 keV. Three model MUSTANG maps were produced using this process and are shown alongside the data in Figure 6. The model uncertainty is dominated by the errors in spectroscopic k_BT_e in region II. To account for this in data comparison, we show three model images corresponding to pressure models produced with the best fit and the temperature fits to Chandra data at ±1σ. The flux scale in the MUSTANG map is completely consistent with the X-ray analysis and is suggestive of a temperature closer to the low end of the allowed 1σ parameter space.

**Figure 6.** MUSTANG SZE data and models in MACS0744. Model images were generated by integrating the best fit and 1σ three-dimensional pressure models from the X-ray along the line of sight and passed through the relevant angular transfer function. Comparison to the MUSTANG data shows the excellent agreement with predicted flux scale. A shock temperature closer to the −1σ value is favored by the measured SZE data.
Download figure:
Standard image High-resolution image

5.1.4. Chandra Spectroscopy

Informed by the shock modeling of the Chandra data, regions corresponding to the core, shock-heated, and pre-shock regions are constructed and used for spectral extraction. These correspond to regions I, II, and III in Figure 4. Since the calibration varies both in time and over the ACIS CCDs, spectra are extracted and response files computed for each of the three observations individually. All three spectra are then fit simultaneously.

XSPEC (Arnaud 1996; Dorman & Arnaud 2001) is used to model the ICM with a Mekal spectrum (Mewe et al. 1985, 1986; Liedahl et al. 1995; Arnaud & Rothenflug 1985; Arnaud & Raymond 1992). In this fit, we account for Galactic extinction and assume the solar abundances of Asplund et al. (2009). The cross sections of Balucinska-Church & McCammon (1992) with an updated He cross section (Yan et al. 1998) are used. The "cstat" statistic, which is similar to the Cash (1979) statistic, is used when modeling the data to properly account for low counts. All three spectra are fit simultaneously to the same plasma model with the abundance fixed to be 0.3 solar in all cases. The normalizations are allowed to float between data sets. The fit is limited to photons within the energy range 0.7–7.0 keV. Best-fit values for the electron temperature and 68% confidence ranges are summarized in Table 3. Though the uncertainty is large in the photon-starved shock-heated region, it is clear that there is a significant increase in temperature in this region compared to the surrounding regions. The uncertainty in the Chandra calibration may affect the temperatures in Table 3, systematically pushing them all in the same direction on order 10%–15% (Reese et al. 2010; Nevalainen et al. 2010). The overall calibration of the X-ray spectra, however, will have a much smaller effect on the temperature jump.

5.1.5. Mach Number

We calculate the Mach number of the potential shock by fitting the Rankine–Hugoniot jump conditions. These relations come from the conditions set by conservation of energy across a boundary. For an extensive review of the physics of shocks in a fluid, see Landau & Lifshitz (1959). The Mach number can be obtained independently by fitting the jump in density from X-ray surface brightness or in temperature as measured by spectroscopy. We use the analytic expressions from Finoguenov et al. (2010) for the Mach number in these two cases:

$\begin{equation} \mathcal {M}_{\rho }=\left[ \frac{ 2\frac{\rho _{2}}{\rho _{1}}}{\gamma +1-(\gamma -1)\frac{\rho _2}{\rho _1}} \right]^{1/2} \end{equation} \tag{ 1 }$

and

$\begin{equation} \mathcal {M} _{T} = \left\lbrace \frac{ 8 \frac{T_2}{T_1} - 7 + \big[ \big(8 \frac{T_2}{T_1} - 7 \big)^2 + 15 \big]^{1/2} }{5}\right\rbrace ^{1/2}, \end{equation} \tag{ 2 }$

where we assume the adiabatic index for a monatomic gas $\gamma =\frac{5}{3}$ and ρ₁, ρ₂, T₁, and T₂ are the density and temperature before and after the shock.

The Mach number can also be calculated from the stagnation condition. This relates the ratio of the pressure at the edge of the cold front, P_st, over the pressure just ahead the shock front, P₁, to the Mach number through the relationship

$\begin{equation} \frac{P_{{\rm st}}}{P_{1}} = \mathcal {M}_{{\rm st}}^{2}\left(\frac{\gamma +1}{2}\right)^{\frac{\gamma +1}{\gamma -1}}\left(\gamma - \frac{\gamma -1}{2\mathcal {M}_{{\rm st}}^{2}}\right)^{-\frac{1}{\gamma -1}} \end{equation} \tag{ 3 }$

as presented in Sarazin (2002).

We calculate the Mach number for the potential merger in MACS0744 using Equations (1), (2), and (3). The value obtained from the density jump conditions was calculated from the posterior MCMC used in the fit to the X-ray surface brightness. This yielded the value $\mathcal {M}_{\rho }=1.2^{+0.2}_{-0.2}$ where the errors are 1σ and the full discrete probability distribution function is shown in Figure 8. The Mach number obtained from the relation imposed by the stagnation condition is $\mathcal {M}_{{\rm st}}=1.4^{+0.2}_{-0.2}$ which is in excellent agreement with the number provided by fitting the density jump. The temperature jump conditions at the shock yield a higher value, $\mathcal {M}_{T}=2.1^{+0.8}_{-0.5}$ . While this measurement suggests a greater shock velocity, the error bars are large and it agrees at the 1.3σ level with our estimate from the density jump condition. The flux scale in the MUSTANG image suggests that the true temperature is toward the low end of the Chandra range as is shown in Figure 6. The shock velocity in this cluster is 1827⁺²⁶⁷_{− 195} km s⁻¹ assuming the Mach number obtained from the density jump conditions.

5.1.6. Discussion

This high-redshift system has proved to be an excellent example of the power of combining resolved SZE and X-ray imaging. The high-resolution SZE measurements reveal a region which is likely the result of a shock. Guided by these data, two sharp discontinuities and a spectrum consistent with a substantially hotter plasma are detected in the low S/N X-ray data. Deeper Chandra observations of this cluster will help confirm the presence of a shock and more accurately determine its Mach number, which for the density jump fit to the current data is mildly consistent with a transonic event ( $\mathcal {M}=1$ ).

Figure 7 shows a composite image of this system including the strong lensing mass distribution (J. Richard et al. 2011, in preparation; see also Jones et al. 2010). This reveals a highly asymmetric elliptical mass distribution elongated to the west consistent with the red sequence member galaxy distribution presented in Kartaltepe et al. (2008). Zitrin et al. (2011) have also done a mass reconstruction and independently obtained a similar mass distribution. The Hubble Space Telescope (HST) data shown in green contain multiple bright red elliptical galaxies with brightest cluster galaxy (BCG)-like characteristics. Galaxy G1 is coincident with the X-ray surface brightness peak and is assumed to be the BCG of the main cluster. Roughly 1 arcmin to the west of the X-ray center, the lensing mass reveals a second peak containing the bright red galaxies G2 and G3. While the baryon distribution is elongated in the direction of this potential subcluster, there is no X-ray peak associated with it. This is likely explained by ram-pressure stripping during passage of the subcluster through the main core. Galaxy G4 is another massive cluster member located west of the main peak. This too has a significant dark matter halo with no apparent baryonic peak. The presence of multiple peaks in dark matter and galaxy density with no accompanying baryonic mass is suggestive of a merger scenario in which a smaller cluster has passed through the main core, stripping it of its baryons and producing a shock wave in the ICM. The geometry of the westerly elongated multiply peaked dark matter distribution is qualitatively suggestive of a merger scenario in which the shock-heated gas identified by MUSTANG could have been produced. However, an accurate interpretation of the merger dynamics requires detailed modeling through hydrodynamical simulations.

**Figure 8.** Mach number in MACS0744 obtained with three methods. The black histogram shows the discrete PDF calculated from the density jump conditions in the MCMC fit to the surface brightness distribution. The red area shows the 68% confidence region from the calculation of a temperature jump across the shock front as measured from *Chandra* spectroscopy. The blue area is the 68% confidence region obtained by fitting the stagnation condition. The solid red and blue lines are the best-fit values obtained from the temperature jump and stagnation conditions, respectively. While the stagnation and density jump conditions yield highly consistent results, the result from the temperature jump conditions appear to be biased high. This is due to a heavy reliance on the spectroscopy in the low S/N region II.
Download figure:
Standard image High-resolution image

5.2. RXJ1347−1145 (z = 0.45)

The rich cluster RXJ1347−1145 is an extremely X-ray luminous galaxy cluster (Schindler et al. 1997; Allen et al. 2002) and has been the object of extensive study in SZE, X-ray, lensing, radio, and optical spectroscopy (e.g., Schindler et al. 1997; Pointecouteau et al. 1999; Komatsu et al. 2001; Allen et al. 2002; Cohen & Kneib 2002; Kitayama et al. 2004; Gitti et al. 2007; Ota et al. 2008; Bradač et al. 2008; Miranda et al. 2008). It is one of only a few systems which have been studied with the SZE on sub-arcminute scales thus far and is an excellent example of the potential for the SZE to reveal the rich substructure exhibited in clusters. Initial measurements made with ROSAT reported by Schindler et al. (1997) deemed RXJ1347−1145 a relaxed system as it showed a round, singly peaked surface brightness morphology and a strong cool core. SZE observations made with the NOBA bolometer camera on the Nobeyama 45 m at 150 GHz (Komatsu et al. 2001; Kitayama et al. 2004) revealed an enhancement to the SZE 20'' (170 kpc) to the southeast of the cluster center. This asymmetry is now supported by X-ray and radio data (Allen et al. 2002; Gitti et al. 2007) and is interpreted as a hot (k_BT_e > 20 keV) feature caused by shock heating in a recent merger event (Kitayama et al. 2004). This feature was confirmed recently by MUSTANG (Mason et al. 2010).

The data we present in this paper are identical to those described in Mason et al. (2010) but are processed with the additional per-pixel high-pass Fourier filter described in Section 3. By deliberately isolating the high-frequency spatial modes, we are able to increase the signal to noise on the small-scale features shown in Figure 9. This comes at the expense of attenuating signals with extents greater than the FOV (see Figure 2), which are not reliably recovered due to the aggressive filtering required to subtract atmospheric noise from our observations. This effective high-pass filtering isolates cluster substructure features on the scales of 8''–42'', which for RXJ1347 corresponds to 50–250 kpc. The MUSTANG maps therefore accentuate shock-heated and adiabatically compressed gas in disturbed clusters, but are not suitable for recovering the average cluster scaling properties typically determined on ∼Mpc scales (see, e.g., Bonamente et al. 2008).

**Figure 9.** Top: MUSTANG map of RXJ1347−1145 at 10'' resolution. Contours are in units of 1σ starting at 3σ. Bottom: signal-to-noise map convolved to 18'' resolution. The peak in the southeast quadrant, originally identified by Nobeyama at 4.2σ (Komatsu et al. 2001), is 13.9σ.
Download figure:
Standard image High-resolution image

The southeast enhancement in RXJ1347 is likely due to shock-heated gas caused by a merger, but the geometry and direction of propagation are not obvious as in the case of MACS0744. This makes it difficult to fit the Rankine–Hugoniot jump conditions across a discontinuity in density inferred from X-ray surface brightness. Komatsu et al. (2001) first reported the southeast enhancement at a peak decrement of 4.2σ. Figure 9 contains a signal-to-noise map of the MUSTANG data set smoothed to comparable resolution to that of the original Nobeyama map. This result is now confirmed at a 13.9σ significance level in SZE at 18'' resolution.

We find good qualitative agreement of our map to a model cluster selected from a suite of hydrodynamical simulations by ZuHone et al. (2010). This particular simulated cluster had recently undergone an off-axis merger with a high mass ratio (∼10: 1). At the epoch of observation, the subcluster is moving through the atmosphere of the main cluster on its first pass of the merger. This model also reproduces several other phenomenological features such as a sharp edge in X-ray surface brightness to the east of the core caused by sloshing of the cold gas, and the location of a second bright elliptical galaxy, thought to be the BCG of the infalling subcluster.

5.3. CL1226.9+3332 (z = 0.89)

With a mass of (1.4 ± 0.2) × 10¹⁵ M_☉ (Jee & Tyson 2009) within r₂₀₀, this system is among the largest known in the high-redshift universe. Early measurements of the baryons were reported by Ebeling et al. (2001b) who identified it in the ROSAT WARPS survey. With the limited resolution of ROSAT, the cluster was deemed to display relaxed morphology.

Because of its high redshift, X-ray spectroscopy on this object is difficult. Initial temperature measurements by Chandra (Bonamente et al. 2006) indicated a hot ICM (∼14 keV). This was consistent with Maughan et al. (2004) who made previous measurements with XMM-Newton (∼12 keV). A more detailed analysis of the ICM properties by Maughan et al. (2007) combined Chandra and XMM-Newton spectroscopy. They confirmed the hot ICM and found an asymmetry in the temperature map with the cluster emission southwest of the cluster center hotter than ambient. This object has also been mapped in the SZE on arcminute scales by the SZA and a strong central decrement was measured (Muchovej et al. 2007; Mroczkowski et al. 2009).

Jee & Tyson (2009) mapped the dark matter distribution of this system through a weak lensing analysis. They found that on large scales the cluster was consistent with a relaxed morphology but the core was resolved into two distinct peaks: the dominant one in close proximity to the BCG and another ∼40'' to the southwest. While this subclump shows no surface brightness peak in either Chandra or XMM-Newton data, the location is consistent with the temperature enhancement reported by Maughan et al. (2007) and a secondary peak in the member galaxy density. One possible explanation of this is a merger scenario in which a smaller cluster has passed through the dominant core on a southwest trajectory stripping its baryons and causing shock heating.

The MUSTANG map is shown in Figure 10. It reveals an asymmetric, multiply peaked pressure morphology in this high-z system. The most pronounced feature is a narrow ridge ∼20'' long located ∼10'' southwest of the X-ray peak. A second peak is found in good proximity to the X-ray emission which is also coincident with the BCG. Also shown in this figure is the X-ray-derived temperature and pseudopressure (defined as the product of the temperature map and the square root of surface brightness) maps from Maughan et al. (2007). The Chandra surface brightness image was produced with 74 ks of archival data taken in ObsIDs 3180, 5014, and 932. There is good qualitative agreement between the two data sets, although the small-scale features seen by MUSTANG are absent in the Maughan map. This discrepancy can be explained by the heavy reliance of the X-ray-derived pseudopressure on the temperature map which was produced with a variable-sized aperture. Therefore, adjacent pixels are not independent. This correlation makes the map less sensitive to small-scale features.

Figure 11 shows radial profiles of the X-ray surface brightness, SZE, and lensing mass distribution. The profiles are centered on the X-ray peak which is coincident with the BCG. Each plot shows two profiles, one taken from the southeastern quadrant (red) and the other in the southwestern quadrant (black). It is clear that all data sets are consistent with an asymmetry elongated toward the southwest as proposed by the merger scenario.

The core of this cluster is compact on the sky due to its high redshift. For this reason, the MUSTANG map can be expected to contain non-negligible amounts of flux that have been modeled in other SZE observations. To quantify the significance of the substructure, we compare our map to the best-fit spherically symmetric Nagai et al. (2007) model of the SZA data as presented by Mroczkowski et al. (2009). Figure 12 shows our map. We assume that the spherically symmetric extended component model is centered on the X-ray peak and takes the difference between the two maps. The residual map figure contains the peak of the ridge at a 4.6σ level.

There is a positive unresolved feature in the residual map in Figure 12 located 8'' northwest of the X-ray peak. This unexplained feature could have several interpretations. It could simply be a noise artifact. However, it is also possible that it is a faint unresolved source. Because it is not detected at 30 GHz (Mroczkowski et al. 2009), such a source would require a rising spectrum in the millimeter as would be expected from a high-redshift, dusty star-forming galaxy. This galaxy could be lensed as speculated by Blain et al. (2002) and Lima et al. (2010) and similar to the one found in the Bullet Cluster by Wilson et al. (2008) and confirmed by Rex et al. (2009). Disentangling speculation such as this requires the addition of resolved millimeter or submillimeter follow-up with different instruments.

A multi-wavelength composite image of this system is presented in Figure 13 which includes the weak lensing mass distribution presented in Jee & Tyson (2009). The northern end of the dominant ridge in the MUSTANG image, labeled "B" in this figure, is roughly coincident with the lensing mass peak. The orientation of the ridge is approximately orthogonal to a vector connecting the BCG and secondary lensing peak which is most likely the trajectory of the subcluster. We posit that the ridge is produced by a reservoir of shock-heated gas created in the core passage of the subcluster, reminiscent of the eastern peak in the famous "Bullet Cluster" (Markevitch et al. 2002).

5.4. A1835 (z = 0.25)

The massive cool core cluster A1835 has proved to be an excellent laboratory for studying a range of cluster physics. It has been used to map the large-scale dark matter distribution (Clowe & Schneider 2002), look for effects of small-scale turbulence (Sanders et al. 2010), search for lensed background submillimeter galaxies (Ivison et al. 2000), map the extended radio emission (Govoni et al. 2009), and study the central cool core (e.g., Peterson et al. 2001; Schmidt et al. 2001). This cluster has also been the subject of extensive SZE modeling (Reese et al. 2002; Benson et al. 2004; Bonamente et al. 2006, 2008; Horner et al. 2010).

Aside from the central ∼10'' region which displays a cavity system excavated by a central AGN (McNamara et al. 2006), the X-ray morphology is well described by a spherically symmetric geometry with no obvious substructure. This distinguishes it from the rest of our sample. The absolute pressure is extremely high in the core, as was demonstrated by Reese et al. (2002) who measured a central decrement of −2.502^+0.150_{− 0.175} mK at 30 GHz. However, the MUSTANG map shown in Figure 14 contains a low signal-to-noise detection of the SZE. The map in this figure has been smoothed to an effective resolution of 18'' to increase the signal to noise. This figure also shows the Chandra image produced from 222 ks of data, merged from ObsIDs 495, 496, 6880, 6881, and 7370. As described in Section 3, the filtering techniques applied are optimized to produce high signal-to-noise maps on small-scale structures. Therefore, the lack of high significance SZE in the reconstructed image is indicative of a featureless, smooth, broad signal.

**Figure 14.** Left: MUSTANG SZE image of A1835 smoothed to 18'' resolution. Contours are units of 0.5σ starting at 2.5σ. Note the central unresolved radio source. Right: *Chandra* 0.7–7.0 keV image of A1835 smoothed with a 15 Gaussian. Contours on the right are identical to those on the left.
Download figure:
Standard image High-resolution image

Figure 15 compares several pixel histograms of the signal region in the MUSTANG A1835 map to asses the significance of any deviation from azimuthal symmetry in the moderate detection. To account for the non-uniform noise distribution caused by the centrally weighted scan patterns used to map the source, each histogram is generated from the product of the map and the square root of the weight map. There are four histograms plotted, each is described below.

1.
The green trace shows the pixel distribution of the circular region centered on the radio point source in the BCG, 1' in diameter. The distribution is far from Gaussian, with the peak significantly below zero, dominated by the SZE signal. While the overall significance of the detection in a given beam is moderate, the region shown is many beams across. The histogram also shows a tail with fewer pixels extending far to the positive. This is dominated by emission from the central radio source.
2.
The black trace is a pixel histogram taken from a region off source which contains negligible signal. This is taken to be indicative of the noise in the map. The 1σ noise level is defined to be the standard deviation of this distribution.
3.
A Gaussian model of the noise is shown as the blue trace. The σ of this Gaussian is equal to the standard deviation of the off-source region shown in black.
4.
An azimuthal average of the map was generated using the location of the radio point source as the centroid. A residual map was then generated by taking the difference of the nominal map and the azimuthally symmetric map. The red trace shows the pixel distribution of the central 1' diameter circular region in the residual map. Obvious deviations from spherical symmetry in the pressure structure of A1835 are expected to reveal a significant departure from a Gaussian distribution in this trace. No significant discrepancies are seen in this distribution. From this we conclude that at the depth of our map no significant pressure substructure is detected.

6. POINT-SOURCE CONTAMINATION

Many clusters, particularly those with cool cores, are known to harbor radio AGNs, often in the central galaxy, with non-negligible luminosities at 90 GHz. These sources contaminate the SZE signal, obscuring the plasma core from a MUSTANG image. In a single dish measurement, the flux of the point source is degenerate with the amplitude of the SZE and therefore difficult to constrain without independent information. In order to account for any compact radio sources in these regions, we rely on archival radio data and extrapolate the flux measurements from FIRST and NVSS at 1.4 GHz and SZA, OVRO, and BIMA at ≈30 GHz (LaRoque et al. 2006; Mroczkowski et al. 2009) to 90 GHz (see Table 4) using a power-law fit, where the uncertainties are treated in quadrature. The predicted fluxes of the radio sources are given in Table 4. Temporal variability in source flux is not accounted for although the multi-band measurements were obtained at different epochs. We are forced to make the assumption of flux stability as there are insufficient available time-sampled radio data to do otherwise. Of the cluster fields presented here, RXJ1347 and A1835 contain significant contributions at 90 GHz from radio point sources clearly detected by MUSTANG. We infer no physics from the central 10'' regions of these clusters as they are obscured by AGN in the BCGs.

Table 4. Unresolved Radio Sources

Cluster Field	Coordinates (J2000)^a		Flux (1.4 GHz)	Flux (30 GHz)^b	Flux (90 GHz)^c
	α	δ	NVSS/FIRST (mJy)	(mJy)	(mJy)
A1835	$14^{\rm h} 01^{\rm m} 02\mbox{$.\!\!^{\mathrm s}$}1$	+02°52'432	31.25 ± 1.57/39.32 ± 1.56	2.8 ± 0.3	1.2 ± 0.2/1.1 ± 0.2
RX J1347.5−1145	$13^{\rm h} 47^{\rm m} 30\mbox{$.\!\!^{\mathrm s}$}7$	−11°45'086	45.89 ± 1.46/NA	10.4 ± 0.3	6.2 ± 0.3
RX J1347.5−1145	$13^{\rm h} 47^{\rm m} 30\mbox{$.\!\!^{\mathrm s}$}1$	−11°45'302	17.66 ± 3.16/NA	<0.3	<0.07
CL J1226.9+3332	$12^{\rm h} 26^{\rm m} 58\mbox{$.\!\!^{\mathrm s}$}2$	+33°32'486	3.61 ± 0.22/4.34 ± 0.47	<0.2	<0.13

Notes. ^aCoordinates are from FIRST except in RXJ1347 where they come from NVSS. ^bMeasured by OVRO, BIMA, and the SZA. ^cExtrapolated assuming a power-law spectral energy distribution.

Download table as: ASCII Typeset image

In the field of RXJ1347 there is a second radio source, ∼1' to the southwest of the central galaxy detected by NVSS as well as Gitti et al. (2007). It is also likely to be associated with a submillimeter source identified by Kitayama et al. (2004) in SCUBA data. This source is not significantly detected by MUSTANG. This is expected as extrapolations from lower frequencies place the predicted flux well below the noise floor. The BCG in CL1226 also contains a radio source measured by NVSS and FIRST. Extrapolations of the spectral index suggest it should be weak at 90 GHz and it is not detected by MUSTANG. In the field of MACS0744, there is no previously reported radio source.

7. SZE FLUX ESTIMATES

In this section, we provide estimates of the SZE flux measured in MUSTANG maps at high significance and on scales sampled by the instantaneous FOV (see Figure 2). We quote the binned flux within two radii, θ_2σ and θ_3σ, which correspond to the radii within which the mean significance is >2σ per beam, and again at >3σ per beam, binned over many beams. We compare our estimates to extrapolations of the SZE flux reported by Bonamente et al. (2008), which was computed within r₂₅₀₀ using 100 kpc core-cut β-model fits to the lower resolution 30 GHz OVRO/BIMA data sets. We use the Itoh et al. (1998) relativistic corrections to the SZE flux frequency relation (described in Appendix B) when scaling the 30 GHz measurements to 90 GHz, assuming for the relativistic correction the isothermal temperatures reported by Bonamente et al. (2008). Any bias due to non-isothermality as well as any discrepancy between temperatures determined with newer calibrations (Reese et al. 2010) and those reported in Bonamente et al. (2008) lead to <2% bias in this frequency rescaling, which is well within the calibration and compact radio source contamination uncertainties in both measurements.

For A1835, CL1226, and RXJ1347, we bin the flux per pixel within a circular region centered on the peak S/N of the map. For MACS0744, we choose an elliptical region to capture, approximately, the shape of the prominent SZE substructure (i.e., the shock-heated region reported in Section 5.1.1). Uncertainties in the absolute calibration are on the order of ≲ 15%, which we include in our estimates. Since we measure SZE flux as a decrement, we add the expected contribution from radio sources (see the flux estimates in Table 4) to our measurements to obtain estimates of the underlying SZE flux reported in Table 5. In the flux estimates presented, we have accounted for attenuation by the filtering by dividing by the mean amplitude of the angular transfer function over the Fourier modes between the beam scale and the radius to which we are integrating. This value is on the order of ∼0.7 from beam (9'') to FOV (42'') scales (see Figure 2).

Table 5. SZE Flux Estimate Comparison

Cluster Name	z_r	D_A	θ₂₅₀₀^a	Y^b	\|F_SZE(90 GHz)\|^c	\|F_{SZE, MUSTANG}\|^d	θ_3σ	θ_2σ	E_th
		(Gpc)	('')	(10⁻¹⁰)	(mJy)	(mJy)	('')	('')	(10⁶² erg)
A1835	0.25	0.81	172 ± ⁵₄	2.09 ± ^0.17_0.16	174 ± ¹⁴₁₃	2.4–3.4±^0.3_0.3	22.9	27.3	0.6–0.9
RX J1347.5−1145	0.45	1.19	122 ± ⁴₄	1.62 ± ^0.18_0.18	135 ± ¹⁵₁₅	12.9–18.5±^2.3_2.3	18.8	22.1	7.1–10.2
MACS J0744.8+3927	0.69	1.47	59 ± ³₃	0.34 ± ^0.04_0.04	28 ± ^3.3_3.3	0.8–1.2±^0.1_0.1	6.8	9.8	0.7–1.0
CL J1226.9+3332	0.89	1.60	66 ± ⁷₆	0.35 ± ^0.05_0.05	29 ± ^4.2_4.2	2.1–2.6±^0.4_0.4	15.0	17.9	2.1–2.6

Notes. ^aThe angular radius on the sky within which the density is 2500 times the critical density of the universe at the cluster's redshift. ^bIntegrated Compton Y = ∫y_C dΩ values from Bonamente et al. (2008). ^cSZE flux measurements on arcminute scales, as reported from fits to OVRO/BIMA data (Bonamente et al. 2008). The 30 GHz, low-resolution measurements were scaled to 90 GHz. ^dFirst and second numbers report the SZE flux binned from pixels within radii where the mean significance is greater than 3σ (θ_3σ) and 2σ (θ_2σ) per beam, respectively.

Download table as: ASCII Typeset image

Measurements of SZE flux provide an estimate, without relying on X-ray data, of the thermal energy in the ICM substructure we see. Table 5 contains the integrated flux and thermal energy estimates for these four objects calculated using the methods described in Appendix B.

It is important to note that the flux and energy estimates reported here from the MUSTANG measurements track pressure substructure on scales ≲ 30''. Contributions from larger scales do not remain after the filtering described in Section 3. We compare these fluxes with those from the smooth, arcminute-resolution SZE signal fitted to the lower resolution OVRO and BIMA 30 GHz SZE data, reported in Bonamente et al. (2008). Assuming that the line-of-sight extents of these structures are comparable to their angular extents, the volumes sampled in the binned regions of the MUSTANG maps are <2% of that within the θ₂₅₀₀ values reported in Bonamente et al. (2008); the measured fluxes from these substructures, however, contribute as much as 10% of the flux within θ₂₅₀₀, implying significant overpressure from within these regions. Since residual emission from the extended signal has not been subtracted, these values should be regarded as upper limits of the energy contained in the small-scale structure.

While A1835 and RXJ1347 have comparable integrated flux on large scales, MUSTANG measures much more signal from the disturbed, merging cluster RXJ1347. The relatively insignificant small-scale (≲ 30'') signal in the A1835 is likely due to the smooth, relaxed nature of this cool-core cluster.

8. CONCLUSIONS

In this paper, we have presented high-resolution images of the SZE in four massive galaxy clusters produced from MUSTANG observations. Three of the four systems probed here display substructure in the SZE. In the case of MACS0744, we identify a likely shock front propagating with a Mach number of $\mathcal {M}=1.2^{+0.2}_{-0.2}$ . The shock-heated kidney-shaped feature is located between the system's main mass peak and a second peak which shows no evidence of significant baryonic mass. In our highest redshift system, CL1226, we find a multiply peaked pressure distribution with an asymmetric morphology. The location and orientation of a ridge found in the SZE, along with a southwesterly elongated shape, are qualitatively supportive of the merger scenario proposed by Jee & Tyson (2009). We also present a new reduction of the data from observations of RXJ1347 presented in Mason et al. (2010). This higher signal-to-noise map confirms the previously reported southeast pressure enhancement at a 13.9σ confidence level. In A1835 we report a detection consistent with a spherically symmetric pressure distribution and no significant substructure.

This pilot study has demonstrated the potential of high-resolution SZE to identify substructures such as weak shocks in galaxy clusters. This is particularly true of the high-redshift universe where the X-ray data are photon starved. A next generation feedhorn-coupled TES bolometer array for the GBT is currently in the planning stages. With a much larger FOV (4 farcm 5) and a mapping speed 1000 times that of MUSTANG, it will be able to image a large number of clusters on angular scales from 9'' to 9'. Other instruments coming online in the next decade, such as ALMA, the LMT, SCUBA2, and CCAT, will also have high-resolution SZE capabilities.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The observations presented here were obtained with telescope time allocated under NRAO proposal IDs AGBT08A056, AGBT09A052, AGBT09C059, and AGBT10A056. We thank James Jee and Anthony Tyson for providing their lensing map for CL1226, Johan Richard and J. P. Kneib for their lensing map in MACS0744, John Zuhone for providing a Compton y_C map from his numerical simulation of RXJ1347, and Ben Maughan for his temperature and pressure maps in CL1226. The contributions of Dominic Benford, Harvey Moseley, Johannes Staguhn, Jay Chervenak, Kent Irwin, Peter Ade, Carole Tucker, Bill Cotton, and Mark Whitehead were essential to the functionality of the instrument. The late night assistance of the GBT operators was much appreciated during the observations. We also thank James Aguirre, Danny Jacobs, and Gary Bernstein for useful conversations, and the anonymous referee for comments that have improved the focus and clarity of the work. Much of the work presented here was supported by NSF grant AST-0607654. P.M.K. was also funded by the NRAO graduate student support program. T.M. is supported as a NASA Einstein Fellow and obtained funding through grant PF0-110077. C.L.S. and M.S. were supported in part by Chandra grants G07-8129X, GO8-9083X, GO9-0135X, and G09-0148X and XMM-Newton grants NNX08AZ34G, NNX08AW83G, and NNX09AQ01G.

APPENDIX A: SHOCK MODEL

A.1. Surface Brightness Profiles

In this work, we measure the density characteristics of a shock front and cold front in MACS0744 by analyzing the elliptical profiles of X-ray surface brightness I(x, y) in some observed photon energy band E₁ to E₂. (In this paper, we consider the surface brightness in the 0.7–7 keV band.) Here, we give the analytic expressions for the X-ray surface brightness of elliptical X-ray images with discontinuities. The X-ray surface brightness is given by the line-of-sight integral

$\begin{equation} I(x,y)= \frac{1}{4\pi (1+z_{r})^{\eta }}\int \varepsilon (x,y,z) \, dz \,, \end{equation} \tag{ A1 }$

where ε is the X-ray emissivity integrated over all directions in the emitted energy band E₁(1 + z_r) to E₂(1 + z_r). The Cartesian coordinates x, y, and z are aligned as shown in Figure 16, with z being along the line of sight. The cluster redshift is z_r. The parameter η = 4 if I(x, y) is given in energy units, and η = 3 if I(x, y) is in counts units, which is generally the case for X-ray observations.

**Figure 16.** Elliptical geometry for a single surface brightness edge used in modeling shock fronts.
Download figure:
Standard image High-resolution image

We fit the data with an analytic expression for the above integral obtained with the following assumptions.

1.
The X-ray emissivity ε(x, y, z) is constant on concentric, aligned, similar ellipsoidal surfaces with the geometry and conventions described in Figure 16. The three principal axes of this elliptical distribution are a, b, and c.
2.
Two of the principal axes of the distribution (a and b) lie in the plane of the sky, and the third axis (c) lies along the line of sight. We take the x-axis of our coordinate system to be parallel to a, and the y-axis to be parallel to b. The axis given by a or x is along the direction of propagation of the shock and/or cold front.
3.
Between each of the discontinuities, the emissivity varies as a power law of the radius, ε = ε₀r^−p. Here, r is the scaled elliptical radius r = [(x/a)² + (y/b)² + (z/c)²]^1/2 and p is the power-law index. The emissivity changes discontinuously at the shock front and/or cold front.
4.
The shock front and/or cold front has rotational symmetry about an axis in the plane of the sky along its direction of propagation (c = b). Although we make this assumption in our analysis of the data on MACS0744, none of the expressions given below depend on this assumption, and are correct for any c.

We treat separately each of the regions bounded by one or two discontinuities. In the case of a shock and cold front, there are three separate regions: the pre-shock gas, the shock-heated gas, and the cold-front gas. Since Equation (A1) is linear in ε, we can then sum the surface brightnesses of these regions to give the total surface brightness.

Since the plane of the sky corresponds to a plane of symmetry at z = 0 in this model, the integral for the surface brightness can be limited to positive z and doubled, giving

$\begin{equation} I(x,y)= \frac{1}{2\pi (1+z_{r})^{\eta }}\int _{q_1}^{q_2} \varepsilon (x,y,z)\,dz \,. \end{equation} \tag{ A2 }$

The values of q₁ ⩾ 0 and q₂ ⩾ 0 give the extent of the cluster region along the line of sight. The general form for the surface brightness for each of the regions obtained with these assumptions after integration is

$\begin{equation} I(x,y)=\frac{1}{4 \pi ^{1/2} (1 + z_r)^\eta } \epsilon _{0}c\frac{\Gamma \big(p-\frac{1}{2}\big)}{\Gamma (p)}A^{-2p+1}\phi \,, \end{equation} \tag{ A3 }$

where we define A to be the two-dimensional scaled elliptical radius

$\begin{equation} A(x,y) \equiv \left(\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}} \right)^{1/2} \,, \end{equation} \tag{ A4 }$

and Γ is the standard Gamma function. The piecewise function ϕ takes a form which depends on the complexity of the model for a given region of interest. Between the discontinuities, each emission region can be treated as having a single outer edge, a single inner edge, or both an inner and an outer edge. For example, in MACS0744, the pre-shock region has a single inner edge, the cold core has a single outer edge, and the shock-heated region has both an inner and an outer edge. The total surface brightness is the sum of these three regions.

For one outer edge, we assume this edge is located at r = 1 in three dimensions and at A = 1 in projection. Then, the bounds on the integral in Equation (A2) are q₁ = 0 and

$\begin{equation} q_2=c \left\lbrace \begin{array}{@{}ll}(1-A^{2})^{1/2}, & A < 1 \\ 0 \, & A \ge 1 \,, \end{array} \right. \end{equation} \tag{ A5 }$

and ϕ takes the form

$\begin{equation} \phi = \left\lbrace \begin{array}{@{}ll}1-I_{A^{2}}\big(p-\frac{1}{2},\frac{1}{2}\big), & A < 1\\ 0, & A \ge 1. \end{array} \right. \end{equation} \tag{ A6 }$

Here, I_x(u, v) is the scaled incomplete beta function I_x(u, v) ≡ B_x(u, v)/B(u, v), B_x(u, v) is the incomplete beta function, and B(u, v) ≡ Γ(u)Γ(v)/Γ(u + v) is the beta function. Note that very efficient algorithms for calculating B(u, v) and I_x(u, v) exist and can be found as intrinsic functions on most computer systems. Alternatively, they are given in Numerical Recipes (Press et al. 1993).

For a single inner edge located at r = R in three dimensions and at A = R in projection, the bounds are

$\begin{equation} q_1 = c \left\lbrace \begin{array}{@{}ll}(R^2-A^{2})^{1/2} \,, & A < R \\ 0 \,, & A \ge R \,, \end{array} \right. \end{equation} \tag{ A7 }$

and q₂ = ∞ and we have

$\begin{equation} \phi = \left\lbrace \begin{array}{@{}ll}I_{\frac{A^2}{R^2}}\big(p-\frac{1}{2},\frac{1}{2}\big)\,, & A < R \\ 1 \,, & A \ge R \,. \end{array} \right. \end{equation} \tag{ A8 }$

It is useful to note that our expression for a single outer edge is mathematically identical to the expression derived in Vikhlinin et al. (2001) but uses the incomplete beta function (which is more convenient numerically) as opposed to the hypergeometric function.

Finally, for a region with two edges, we will take their locations to be r = 1 in three dimensions and A = 1 in projection for the inner edge, and r = R or A = R for the outer edge, where R > 1. The bounds on the integral become

$\begin{equation} q_1 = c \left\lbrace \begin{array}{@{}ll}(1-A^{2})^{1/2} \,, & A < 1 \\ 0 \,, & A \ge 1 \, \end{array} \right. \end{equation} \tag{ A9 }$

and

$\begin{equation} q_2=c \left\lbrace \begin{array}{@{}ll}(R^2-A^{2})^{1/2} \,, & A < R \\ 0 \,, & A \ge R \,. \end{array} \right. \end{equation} \tag{ A10 }$

The expression for ϕ becomes

$\begin{equation} \phi = \left\lbrace \begin{array}{@{}ll}I_{A^{2}}\big(p-\frac{1}{2},\frac{1}{2}\big) - I_{\frac{A^{2}}{R^{2}}}\big(p-\frac{1}{2},\frac{1}{2}\big) \,, & A < 1 \\ 1-I_{\frac{A^{2}}{R^{2}}}\big(p-\frac{1}{2},\frac{1}{2}\big) \,, & 1 \le A < R \\ 0 \,, & A \ge R \,. \end{array} \right. \end{equation} \tag{ A11 }$

A.2. Density Profiles

Once we have obtained the power-law index p and the normalization ε₀ by fitting Equation (A3) to the data, we can reconstruct the intrinsic emissivity distribution. This is related to the density distribution n_e(r) by

$\begin{equation} n_{e}(r)=\left[ \frac{\varepsilon (r)}{\Lambda (T_e,Z)} \right]^{1/2} \,, \end{equation} \tag{ A12 }$

where Λ is the X-ray emissivity function which depends on electron temperature T_e and abundance Z.

If XSPEC⁷ is used to determine the temperatures in the emission regions, the same models can easily be used to determine the value of Λ. This has the great advantage that the models, temperature, abundances, and instrument responses used for the spectral analysis will be completely consistent with those used to determine n_e(r). We assume here that the model is a single-temperature MEKAL or APEC model. For this purpose, only the shape of the spectrum matters, not its normalization, so the region fit in XSPEC need not be identical to the region fit in the surface brightness analysis, as long as the spectral shape is assumed to be the same. If the surface brightness I is analyzed in energy units, then the procedure is to determine the X-ray flux F of the spectral region in the same band and with the same instrument as used to fit the surface brightness. If the surface brightness was corrected for absorption, then the absorbing column should first be set to zero. One also needs to record the normalization of the thermal model, which is defined as

$\begin{equation} K \equiv \frac{10^{-14}}{4 \pi (1 + z_r)^2 D_A^2} \int n_e n_p \, dV \,, \end{equation} \tag{ A13 }$

where D_A is the angular diameter distance to the cluster, n_p is the proton number density, and V is the volume of the emitting region. Then, the relevant X-ray emissivity function is

$\begin{equation} \Lambda = \frac{F(1 + z_r)^2 }{10^{14}\, {\rm K} (n_e / n_p) } \,. \end{equation} \tag{ A14 }$

Here, n_e/n_p ≈ 1.21 is the ratio of the electron-to-proton number densities, and is essentially a constant for typical cluster temperatures and abundances.

If the surface brightness is determined in count units (as is typically the case with X-ray observations), then the procedure is to set the observed energy band and instrument in XSPEC to the one used for the surface brightness measurements, and then type "show" to determine the model count rate, CR. Then, the emissivity function is

$\begin{equation} \Lambda = \frac{{\rm CR} \,(1 + z_r)}{10^{14}\, {\rm K} (n_e / n_p) } \,. \end{equation} \tag{ A15 }$

A.3. Pressure and SZE

With a three-dimensional density model obtained through the above procedure and measurements of T_e from X-ray spectroscopy, one can produce a three-dimensional pressure model which can be used to predict the observed SZ flux. From the ideal gas law, the electron pressure is simply

$\begin{equation} P_{e}(r)=k_{B}n_{e}(r) T_e (r) \,. \end{equation} \tag{ A16 }$

By integrating this along the line of sight, one can obtain a two-dimensional map of the Compton y_C parameter

$\begin{equation} y_C (x,y)=\int \frac{P_{e}(r)\sigma _{T}}{m_{e}c_{l}^{2}}\,dz. \end{equation} \tag{ A17 }$

Here, k_B is Boltzmann's constant, σ_T and m_e are the Thomson cross section and mass of the electron, respectively, and c_l is the speed of light.

Assuming that T_e(r) is either a constant or is a power-law function of the radius within each region, the electron pressure will vary as a power law of the elliptical radius, and the same analytic expression (Equation (A3)) can be used to determine y_C(x, y). One simply makes the substitution

$\begin{equation} \frac{1}{4\pi (1+z_{r})^{\eta }} \varepsilon \rightarrow \frac{P_{e}(r)\sigma _{T}}{m_{e}c_{l}^{2}} \,. \end{equation} \tag{ A18 }$

From a map of y_C, it is straightforward to produce a model SZE image.

APPENDIX B: THERMAL ENERGY FROM SZE

The surface brightness of a cluster due to the thermal SZE can be expressed, for dimensionless frequency $x_{\nu } \equiv h\nu /\mbox{$k_{{\rm B}}$}\mbox{$T_{{\rm CMB}}$}$ where h is Planck's constant, ν is the frequency, and T_CMB is the primary CMB temperature, as the change ΔI_SZE relative to the primary CMB surface brightness normalization I₀, as

$\begin{eqnarray} \frac{\Delta I_{\rm SZE}}{I_0} &=& \frac{\mbox{$\sigma _{{\rm T}}$}}{m_e c_l^2} \int \!\! g(x_{\nu },T_e) \, \mbox{$k_{{\rm B}}$}n_e T_e \,dz \end{eqnarray} \tag{ B1 }$

$\begin{eqnarray} \hspace{6pt}&=& \frac{\mbox{$\sigma _{{\rm T}}$}}{m_e c_l^2} \int \!\! g(x_{\nu },T_e) \, P_e \,dz \end{eqnarray} \tag{ B2 }$

$\begin{eqnarray} \hspace{-3pc}&\equiv & g(x_{\nu },T_e) y_C. \end{eqnarray} \tag{ B3 }$

The primary CMB surface brightness normalization (in units of flux per solid angle) is $I_0 = 2 (\mbox{$k_{{\rm B}}$}\mbox{$T_{{\rm CMB}}$})^3(h c_l)^{-2} = 2.7033 \times 10^8 \,\rm Jy \,sr^{-1}$ (see, e.g., Carlstrom et al. 2002). The factor g(x_ν, T_e) encapsulates the SZE flux spectral dependence, which is a function of electron temperature when relativistic corrections are taken into consideration. In the classical physics limit,

$\begin{equation} g(x_{\nu }) = \frac{x_{\nu }^4 e^{x_{\nu }}}{(e^{x_{\nu }}-1)^2} \left(x_{\nu } \frac{e^{x_{\nu }} + 1}{e^{x_{\nu }} - 1} - 4\right). \end{equation} \tag{ B4 }$

We integrate the SZE surface brightness in Equation (B1) to relate the SZE flux from a region of the sky to the underlying electron pressure in the measured ICM feature. The integrated SZE flux is computed (using Equations (B2), (B3), and (A17)) as

$\begin{eqnarray} F_{\rm SZE} &=& \int \! I_{\rm SZE} \, d\Omega = g(x_{\nu }) {I_0} \int \! y_C \, d\Omega\nonumber\\ &=& \frac{\mbox{$\sigma _{{\rm T}}$}}{m_e c_l^2} \, g(x_{\nu }) {I_0} \int \! d\Omega \int \! P_e \,dz. \end{eqnarray} \tag{ B5 }$

Since dΩ = dℵ/D²_A(z_r), where dℵ is the area integration element, F_SZE physically relates to the thermal energy E_th content of the gas within a cylindrical region of a cluster (of volume ΔℵΔz). The electron pressure P_e relates to the total pressure P_tot by the electron weighting factor μ_e ≈ 1.17 (assuming standard abundances) as P_tot = (1 + 1/μ_e)P_e. In terms of the flux (Equation (B5)), the thermal energy content is

$\begin{equation} E_{\rm th} =\frac{3}{2} \frac{(1 + 1/\mu _e) \, m_e c_l^2 \, F_{\rm SZE} \, D_A^2(z_r)}{\mbox{$\sigma _{{\rm T}}$}I_0 \, g(x_{\nu })}. \end{equation} \tag{ B6 }$

For MUSTANG data at 90 GHz, and an assumed $\mbox{$k_{{\rm B}}$}T_e = 10 \,\rm keV$ and μ_e = 1.17, this is

$\begin{equation} E_{\rm th} = |F_{\rm SZE}| D_A^2(z_r) \rm \left[\frac{3.9\times 10^{55} \,erg}{mJy \, Mpc^2}\right] \end{equation} \tag{ B7 }$

for D_A(z_r) in Mpc. In this work, we use the Itoh et al. (1998) relativistic corrections to Equation (B4).

MUSTANG HIGH ANGULAR RESOLUTION SUNYAEV–ZEL'DOVICH EFFECT IMAGING OF SUBSTRUCTURE IN FOUR GALAXY CLUSTERS

Article metrics

Permissions

Author e-mails

Author affiliations

Author notes

Dates

ABSTRACT

1. INTRODUCTION

2. INSTRUMENT AND OBSERVATIONS

2.1. MUSTANG

2.2. Observations

3. MUSTANG DATA REDUCTION

4. CHANDRA DATA REDUCTION

5. RESULTS

5.1. MACS J0744.8+3927 (z = 0.69)

5.1.1. MUSTANG Data

5.1.2. Chandra Data

5.1.3. X-Ray Surface Brightness Shock Modeling

5.1.4. Chandra Spectroscopy

5.1.5. Mach Number

5.1.6. Discussion

5.2. RXJ1347−1145 (z = 0.45)

5.3. CL1226.9+3332 (z = 0.89)

5.4. A1835 (z = 0.25)

6. POINT-SOURCE CONTAMINATION

7. SZE FLUX ESTIMATES

8. CONCLUSIONS

APPENDIX A: SHOCK MODEL

A.1. Surface Brightness Profiles

A.2. Density Profiles

A.3. Pressure and SZE

APPENDIX B: THERMAL ENERGY FROM SZE

Footnotes

MUSTANG HIGH ANGULAR RESOLUTION SUNYAEV–ZEL'DOVICH EFFECT IMAGING OF SUBSTRUCTURE IN FOUR GALAXY CLUSTERS

Article metrics

Permissions

Share this article

Author e-mails

Author affiliations

Author notes

Dates

ABSTRACT

1. INTRODUCTION

2. INSTRUMENT AND OBSERVATIONS

2.1. MUSTANG

2.2. Observations

3. MUSTANG DATA REDUCTION

4. CHANDRA DATA REDUCTION

5. RESULTS

5.1. MACS J0744.8+3927 (z = 0.69)

5.1.1. MUSTANG Data

5.1.2. Chandra Data

5.1.3. X-Ray Surface Brightness Shock Modeling

5.1.4. Chandra Spectroscopy

5.1.5. Mach Number

5.1.6. Discussion

5.2. RXJ1347−1145 (z = 0.45)

5.3. CL1226.9+3332 (z = 0.89)

5.4. A1835 (z = 0.25)

6. POINT-SOURCE CONTAMINATION

7. SZE FLUX ESTIMATES

8. CONCLUSIONS

APPENDIX A: SHOCK MODEL

A.1. Surface Brightness Profiles

A.2. Density Profiles

A.3. Pressure and SZE

APPENDIX B: THERMAL ENERGY FROM SZE

Footnotes