The Dust and Molecular Gas in the Brightest Cluster Galaxy in MACS 1931.8-2635

Our ALMA Band 3 observations were obtained in cycle 4 (Project ID 2016.1.00784.S). Our ALMA Bands 6 and 7 observations, including Atacama Compact Array (ACA) observations, were obtained in cycle 5 (Project ID 2017.1.01205.S). Basic information about these observations are summarized in Table 1.

The Band 3 observations combined two ALMA 12 m array configurations, which yielded a 15–3700 m baseline range and a 13,366 s total integration time. The minimum possible synthesized beam and the recoverable angular scale were 054 and 188, respectively. We configured our Band 3 spectral windows (SPWs) to center on sky frequencies of 85.3 GHz, 87.0 GHz, 97.0 GHz, and 99.0 GHz. The spectral channel width was 3.9 MHz (i.e., spectral resolution was 7.8 MHz), which corresponds to ∼11.7 km s⁻¹ velocity width.

We carried out one set of ALMA 12 m array observations at Band 6, with a baseline range of 15–500 m. It was complemented with a set of ACA observations with an identical spectral setup. In these observations, the Band 6 SPWs were configured to be centered on sky frequencies of 239.5 GHz, 241.5 GHz, 255.5 GHz, and 257.5 GHz. The spectral channel width was 15.625 MHz (∼18.1 km s⁻¹). The on-source times for the 12 m array and the ACA observations were 726 s and 1724 s, respectively. Combining these observations yielded a 074 minimum possible synthesized beam and a 204 largest recoverable angular scale.

Similarly, we obtained one set of ALMA 12 m array observations at Band 7, with a baseline range of 15–314 m, paired with the complementary ACA observations. We configured the Band 7 SPWs to center on sky frequencies of 339.5 GHz, 341.4 GHz, 351.5 GHz, and 353.5 GHz. The spectral channel width was 15.625 MHz (∼13.2 km s⁻¹). The on-source times for the 12 m array and ACA observations were 2722 s and 6804 s, respectively. Combining these observations yielded a 081 minimum possible synthesized beam and a 146 largest recoverable angular scale.

We use the 16 bands of Hubble Space Telescope (HST) photometry for the MACS J1931.8-2635 BCG that were obtained between rest wavelengths of 1745 Å and 1.14 μm as part of the CLASH program. These observations are detailed in Postman et al. (2012) and the reduced data set is described in Fogarty et al. (2015, 2017).

Cycle 4 observations were reduced using the CASA software package version 4.7.2 (pipeline r39732); cycle 5 observations were reduced using CASA version 5.1.1–5 (pipeline version r40896) (McMullin et al. 2007). After regenerating the pipeline-calibrated measurement sets, we used the task concat to concatenate the calibrated measurement sets at each band. We used the task plotms to visually inspect the spectra of our target source, and then used the task uvcontsub to generate the continuum data from the line-free spectral channels and to produce continuum-subtracted CO line data.

We independently created images with the Band 3, 6, and 7 line-free continuum data using the multi-frequency synthesis mode of the clean task. The visibilities were naturally weighted, and the image pixel size was 0075. We produced clean masks using clean in interactive mode to identify regions with significant source emission.

Flux in the Band 3 continuum is primarily due to the AGN. Flux in the Band 6 and 7 continuum images is due to thermal emission from extended dust and synchrotron and hot dust emission from the AGN. The AGN in our ALMA data is point-like and is most prominent in Band 3. To yield the dust continuum images, we subtracted the point source in each band using the image arithmetic task fitcomponents to fit a Gaussian model to the point source with a uniform background level. The full widths at half maximum along the major and minor axes and the position angle of the Gaussian were fixed based on the synthesized beam of each image. We fit the peak coordinates of the point source, along with the source flux and background. The coordinates of the point source were the same in both bands. The point source flux at Band 6 was 4.19 ± 0.42 mJy, roughly equal to the 4.28 ± 0.43 mJy measured in Band 7, ruling out the possibility that the emission is primarily due to graybody emission by hot dust in the AGN torus or unresolved circumnuclear region.

We ran clean in "velocity" mode on each continuum-subtracted line data set, using natural weighting, a 0075 pixel size, and a velocity channel width of 25 km s⁻¹. The achieved synthesized beams for Band 3, Band 6, and Band 7 were 082 × 053, 087 × 072, and 094 × 078. In cases where we needed to produce data cubes with identical angular resolution from different data sets, we selected a common synthesized beam that encompassed the natural restoring beam for each band and smoothed using imsmooth with targetres = True. To study the extended, faint molecular gas structures around MACS 1931 BCG, we also produced uv-tapered CO image cubes to yield better brightness temperature sensitivities.

All images or image cubes used in our scientific analysis were primary beam corrected.

Our analyses of the HST and ALMA data require that the images be aligned to a common astrometric reference frame. The astrometric reference frames for both HST and ALMA are calibrated to the International Celestial Reference System. To assess any systematic offsets in the sky positions of features in our field, we cross-matched our HST astrometry to entries from the Gaia DR2 catalog (Gaia Collaboration et al. 2018) and found a mean offset of 95 mas with an rms of 26 mas. The maximum astrometric error in our ALMA data is estimated to be approximately 5% of the synthesized beamwidth, which corresponds to 28 mas, 37 mas, and 41 mas, respectively, for bands 3, 6, and 7. For our study of the BCG in MACS 1931, any cross-comparisons of aligned features in our HST and ALMA data are performed on angular scales larger than 200 mas. This ensures that any systematic astrometric alignment errors between the HST and ALMA data sets will not be a significant concern.

Our ALMA campaign obtained spatially resolved line emission from the J = 1 to 0 transition of CO(1−0) as well as CO(3−2) and partial coverage of CO(4−3). The spatially integrated line spectra and the velocity-integrated intensity maps are shown in Figures 1 and 2, respectively. We measured spectra in regions where the total intensity is detected at the 2σ level or greater. We also re-ran imaging for CO(1−0) and CO(3−2) using a 15 uv-taper, which down-weighted noisier, higher spatial frequency baselines. Using the original and the uv-tapered images enabled us to resolve the compact and luminous knots of CO emission in the core of the BCG and the extended and faint emission features along the Hα "tail" extending to the northwest. Total fluxes in the uv-tapered images are consistent with the un-tapered images.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Spatially integrated flux from the emission lines were fit to Lorentzian profiles, shown as the solid blue lines in Figure 1. All three lines have two velocity components that are redshifted at ∼0 km s⁻¹ and ∼300 km s⁻¹ relative to the optical emission lines in the MACS 1931 BCG.¹⁰ As shown below, the more redshifted component is due in part to a fraction of the molecular gas flowing either inward or outward from the galaxy center, and in part to motions of gas clumps in the center of the galaxy. Based on the best model fit, the CO(1−0) line has an integrated flux of 3.45 ± 0.47 Jy km s⁻¹, the CO(3−2) line has a flux of 29.0 ± 2.8 Jy km s⁻¹, and the CO(4−3) line has a flux of 33.7 ± 3.5 Jy km s⁻¹. Best-fit parameters for the lines are summarized in Table 2.

Compact knots in the naturally weighted CO(1−0) intensity image trace UV knots present in HST WFC3 F225W-F390W photometry (see Figure 2). These knots are also visible in the naturally weighted CO(3−2) image but are less apparent given the shorter baselines and consequently lower spatial resolution of the image. The faint, extended component of the CO(1−0) and CO(3−2) emission traces the distribution of Hα–N ii emission inferred by the broadband difference image from Fogarty et al. (2015). This faint component forms both a diffuse nebula around the dense knots and a "tail" extending ∼30 kpc to the northwest of the optical brightness peak of the BCG (hereafter referred to as the BCG core). To within the angular resolution of our CO imaging, the HST and ALMA data demonstrate a spatial correlation between the atomic and molecular phases of gas in the BCG, and a correlation between dense knots of molecular gas and sites of star formation, suggesting star formation in the BCG is fueled by multiphase material.

The velocity structure of the molecular gas is complex. We observe a systematic offset of ∼300 km s⁻¹ between the gas in the tail and the gas in the BCG core. This offset is clearly seen in the velocity map of CO(3−2) (see Figure 3) and is consistent with the relative velocities between the BCG core and tail obtained by fitting both CO(1−0) and CO(3−2) in the uv-tapered data cubes. In the MACS 1931 BCG the molecular gas shows no clear morphological symmetry around the core or connection to X-ray cavities that would suggest jet-driven molecular gas outflows like those seen in A1835 and A1664 (McNamara et al. 2014; Russell et al. 2014). Hlavacek-Larrondo et al. (2012) detect probable cavities in the core of MACS 1931; these detected features are found to the east and west of the BCG, as opposed to the approximately north–south orientation of the extended molecular gas.

Zoom In Zoom Out Reset image size

Molecular gas in the core is concentrated in bright knots (which we have divided up into regions A–E in Figure 4) with CO(1−0) velocities ranging from −98 ± 11 to 258 ± 18 km s⁻¹ (see Table 3). One of these knots has two velocity peaks, at −49 ± 14 and 258 ± 18 km s⁻¹, so it is either rotating or consists of two knots projected on top of each other along the line of sight. Spectra for these knots are shown in Figure 4, along with the best-fit Lorentzian profiles. The CO(3−2) and CO(1−0) velocities are consistent on the positions of the knots (see Figure 3).

Zoom In Zoom Out Reset image size

The velocity dispersion maps of CO(1−0) and CO(3−2), estimated from analysis of the second moments of the emission lines, are shown in Figure 3. Velocity dispersion distributions for the two lines are similar, and range between ∼25 and ∼400 km s⁻¹. The range of velocity dispersions we observe is similar to observations of CO in other BCGs (Russell et al. 2014, 2016; Vantyghem et al. 2016, 2018). Figure 5 shows the distribution of the observed CO(3−2) velocity dispersions alongside predictions for simulated BCGs with conditions similar to MACS 1931 when the simulated data are smoothed to spatial resolutions approaching that of our ALMA data. The physical beam sizes of our CO images are 4.1 × 2.6 kpc for the CO(1−0) image and 4.3 × 3.6 kpc for the CO (3−2) image. As the beam size used to smooth the simulated data increases and approaches the beam size of the actual observed data, the typical velocity dispersion measured in the simulated data set comes into agreement with our observations.

Zoom In Zoom Out Reset image size

Individual giant molecular clouds in the Milky Way have typical sizes of 5 to 200 pc (Murray 2011). Therefore, the angular resolution of our ALMA observations is insufficient to probe the internal velocity dispersion of individual molecular clouds. Instead, the velocity dispersions we measure are likely dominated by random bulk motion of different clouds that are spatially close to one another on scales smaller than our beam size. In the nearby Virgo cluster, the ALMA-measured velocity dispersion of the molecular cloud in M87 is 27 ± 3 km s⁻¹ (Simionescu et al. 2018). This is likely still under-resolved as the velocity dispersion of resolved molecular clouds in the Milky Way and nearby star-forming galaxies is typically less than 10 km s⁻¹ (e.g., Stark 1984; Brouillet et al. 1998).

Figures 3 and 5 nonetheless show evidence for structure in the velocity dispersion of CO. Dispersions are lowest near the center of the core (i.e., near the position of knot "D" in Figure 4), and increase toward both the southern and northern peripheries of the molecular gas nebula. Inspection of the 1D spectra at several locations with high-velocity dispersions shows evidence of multiple velocity peaks, suggesting that the increased dispersion is almost certainly due to multiple velocity components along the line of sight, possibly as a result of the molecular clouds having a lower filling factor near the outskirts of the CO emission.

We estimate the total molecular gas mass of MACS 1931 to be (1.9 ± 0.3) × 10¹⁰ M_⊙ based on the CO(1−0) intensity in Jy km s⁻¹, I_CO. We calculated the gas mass using the relation

$$\begin{eqnarray}&&{M}_{\mathrm{mol}}=1.05\times {10}^{4}\displaystyle \frac{{X}_{\mathrm{CO}}}{2\times {10}^{20}\tfrac{{\mathrm{cm}}^{-2}}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}}}\displaystyle \frac{{D}_{L}^{2}}{1+z}{I}_{\mathrm{CO}}\end{eqnarray} \tag{ 1 }$$

where D_L is the luminosity distance in Mpc and X_CO is the CO-to-H₂ conversion factor in units of $\tfrac{{\mathrm{cm}}^{-2}}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}}$ (Bolatto et al. 2013). We adopted the value ${X}_{\mathrm{CO}}=0.4\times {10}^{20}\tfrac{{\mathrm{cm}}^{-2}}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}}$ for starbursts and ultra-luminous infrared galaxies (ULIRGs) used in Russell et al. (2017a) to estimate the molecular gas mass in the Phoenix cluster, in order to be consistent when comparing the molecular gas masses of BCGs (by adopting the Galactic value of X_CO, as was done in McNamara et al. 2014, our estimate of the molecular gas mass rises to (9.4 ± 1.3) × 10¹⁰ M_⊙). The gas reservoir we observe in MACS 1931 is approximately equal to the reservoir in the Phoenix cluster ((2.1 ± 0.3) × 10¹⁰ M_⊙, Russell et al. 2017a), which along with a handful of other recent observations (A1835: McNamara et al. 2014; RXJ0821+0752: Vantyghem et al. 2017, 2019; RXC J1504-0248: Vantyghem et al. 2018), make it among the largest molecular gas masses observed in a cool-core BCG.

Images of the dust continuum emission in our Band 6 and 7 observations are shown in Figure 6. These images are shown after performing the AGN point source subtraction described in Section 3.1. Extended dust emission is clearly visible in both bands in the core and tail regions described in Section 4.1. The Band 6 continuum is centered on 336.0 GHz (892 μm) in the rest frame and the Band 7 on 468.5 GHz (640 μm), so the emission we observe in these bands is either due to very cold dust or is the Rayleigh–Jeans tail of the graybody emission of hotter dust. The dust emission in MACS 1931 is most prominent in the BCG core and in the northern tip of the tail, corresponding to the edge of the extended Hα and UV emission seen in the CLASH UV–optical images shown in Figure 2.

Zoom In Zoom Out Reset image size

We checked to ensure that the emission features we observe are not due to free–free emission or a faint, extended synchrotron source. We examined the Band 3 continuum image after convolving it with the Band 6 synthesized beam, since both free–free and synchrotron emission will appear brighter in Band 3. The Band 3 fluxes have detection significances of ≲0.5σ in the apertures labeled Regions 1, 2, and 3 in Figure 6. The 3σ limit in Band 3 is 0.023 mJy/beam, which is about a factor of 3 fainter than the limits in either Bands 6 or 7, implying that any contamination by extended free–free or synchrotron emission in our dust images is below the rms noise in Bands 6 and 7.

We also estimated the total free–free and synchrotron emission in Band 6 due to star formation using the relationship derived in Schmitt et al. (2006). An SFR of 250 M_⊙ yr⁻¹ yields a total flux in Band 6 of 0.08 mJy. If this flux is evenly distributed throughout the region where Band 6 emission is detected at ≥3σ, we predict a Band 6 flux due to star formation of 0.01 mJy/beam, which is well below the rms noise of 0.020 mJy/beam. The predicted Band 3 flux, meanwhile, is 0.014 mJy/beam, which is consistent with being below the detected limit in our Band 3 image. Since the dust emission is higher in Band 7 than Band 6, while the free–free and synchrotron component is lower, we did not perform the same check for Band 7.

Since we have multiple bands of ALMA far-IR photometry, we are able to constrain the dust temperature using a simple model for the flux. Given the long rest wavelengths of our Bands 6 and 7 photometry, we can only reliably derive the temperature of very cold T_d ≲ 30 K dust from the ALMA data alone, since emission from hotter dust in these bands is proportional to both the dust luminosity and T_d. However, we can still place constraints on the temperature of the cold dust in the tail, especially in regions with relatively high flux in Band 6.

Casey (2012) describes a mid-IR power-law + graybody model that represents the non-polycyclic aromatic hydrocarbon continuum emission due to dust in a ULIRG-like environment with the expression

$$\begin{eqnarray}&&S\left(\lambda \right)=N\displaystyle \frac{{\left(c/\lambda \right)}^{3}\left(1-{e}^{-{\left({\lambda }_{0}/\lambda \right)}^{\beta }}\right)}{{e}^{{hc}/\lambda {{kT}}_{d}}-1}+{N}_{{pl}}{\lambda }^{\alpha }{e}^{-{\left(\lambda /{\lambda }_{c}\right)}^{2}},\end{eqnarray} \tag{ 2 }$$

where $S\left(\lambda \right)$ is the flux density in Jy, λ_c is the turnover wavelength separating the regimes where the mid-IR power-law and the graybody dominate the emission, λ₀ ≡ 200 μm, N is the overall normalization, T_d is the dust temperature, α is the mid-IR power-law index, and β is the graybody emissivity index. Both λ_c and the power-law normalization, N_pl, are functions of the other parameters in the model (Casey 2012). Given the long wavelengths of our two ALMA bands, we choose to neglect the mid-IR power-law component of the model in Equation (2) and model the dust emission as a simple graybody.

Since the dust components we will be able to constrain with this method are cold, we need to account for the impact of the cosmic microwave background (CMB) on our observations, which can potentially bias our estimated temperature constraints by several K. We note, however, that this correction does not qualitatively change the results of our analysis. We modify the first term in Equation (2) to be

$$\begin{eqnarray}\begin{array}{rcl}S\left(\lambda \right) & = & N\left[\displaystyle \frac{{\left(c/\lambda \right)}^{3}\left(1-{e}^{-{\left({\lambda }_{0}/\lambda \right)}^{\beta }}\right)}{{e}^{{hc}/\lambda {{kT}}_{d}}-1}\right.\\ & & +\ \left.\displaystyle \frac{{\left(c/\lambda \right)}^{3}{e}^{-{\left({\lambda }_{0}/\lambda \right)}^{\beta }}}{{e}^{{hc}/\lambda {{kT}}_{\mathrm{CMB}}}-1}\right],\end{array}\end{eqnarray} \tag{ 3 }$$

where ${T}_{\mathrm{CMB}}=2.73\left(1+z\right)=3.69$ K is the temperature of the CMB at z = 0.3525. This equation requires one further modification since the interferometer configuration we use does not have baselines short enough to detect the spatial scales that contain the CMB power. Therefore, while we observe the dust absorption of the CMB, we do not observe the CMB itself, and need to subtract an additional term of ${\left(c/\lambda \right)}^{3}{\left({e}^{{hc}/\lambda {{kT}}_{\mathrm{CMB}}}-1\right)}^{-1}$ from Equation (3) (see, e.g., Zhang et al. 2016). The final expression for the Band 7 to Band 6 continuum flux ratio is therefore

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{S\left({\lambda }_{7}\right)}{S\left({\lambda }_{6}\right)} & = & {\left(\displaystyle \frac{{\lambda }_{6}}{{\lambda }_{7}}\right)}^{3}\left(\displaystyle \frac{1-{e}^{-{\left({\lambda }_{0}/{\lambda }_{7}\right)}^{\beta }}}{1-{e}^{-{\left({\lambda }_{0}/{\lambda }_{6}\right)}^{\beta }}}\right)\\ & & \times \ \left(\displaystyle \frac{\tfrac{1}{{e}^{{hc}/{\lambda }_{7}{{kT}}_{d}}-1}-\tfrac{1}{{e}^{{hc}/{\lambda }_{7}{{kT}}_{\mathrm{CMB}}}-1}}{\tfrac{1}{{e}^{{hc}/{\lambda }_{6}{{kT}}_{d}}-1}-\tfrac{1}{{e}^{{hc}/{\lambda }_{6}{{kT}}_{\mathrm{CMB}}}-1}}\right),\end{array}\end{eqnarray} \tag{ 4 }$$

where λ₇ is the rest wavelength of the Band 7 continuum (640 μm) and λ₆ the rest wavelength of the Band 6 continuum (892 μm).

The contour plot in Figure 7 shows the predicted values of the flux ratio as a function of T_d and β. The parameter β depends on the impact of the physical properties of the dust on its emissivity and for local sources has a value between 0.8 and 2.4 (Dupac et al. 2003). Casey (2012) measured β = 1.6 ± 0.38 for a sample of ULIRGs, which we take to be the most plausible range for β in the MACS 1931 BCG.

Zoom In Zoom Out Reset image size

We examined three regions outside the core that are relatively bright in Band 6, denoted by the solid ellipses in Figure 6 and labeled Regions 1–3. We smoothed the Band 6 and 7 images to a common beam, and measured the flux ratios in these regions. The Band 7 to Band 6 flux ratios for these regions are 1.25 ± 0.53, 3.29 ± 0.99, and 1.29 ± 0.33, respectively. We also measure the flux ratio in the BCG core, obtaining a value of 2.74 ± 0.46, and the overall ratio in the tail, finding a value of 2.00 ± 0.40.

Shaded bands showing the 1σ intervals around the ratios for Regions 1–3 are shown in Figure 7. There is a very cold dust component in Regions 1 and 3; for 1.22 ≤ β ≤ 1.98, we estimate T_d ≲ 8.6 K, and ≲7.3 K at the 1σ confidence level, and ≲17.3 K, and ≲10.2 K at the 95% confidence level respectively. If β = 1, the temperatures for these regions are ≲26.4 K and ≲12.5 K at 95% confidence.

A lower limit of 8.4 K can be placed on the dust temperature of Region 2. If we assume β = 1.6, then we can constrain T_d to >11.2 K. The measured value of the flux ratio of 3.29 in Region 2 implies a temperature of 30.5 K if β = 1.98, but >70 K once β < 1.8, so depending on the emissivity the dust in this region may either be relatively cold or warm.

In the core of the BCG, if we assume β = 1.6, we can place a constraint of T_d > 11.0 K given the 1σ error bars on the flux ratio. However, the tail seems to be substantially colder overall; we find that 5.0 < T_d < 21.9 K for β = 1.6 ± 0.38.

We note that when we model and subtract the bright point sources in the Band 6 and 7 continuum images, the peak positions of the best-fit point source in the two bands differ by less than 0.25 pixels, so the low flux ratios we obtain are not likely to be the result of astrometric offsets between bands.

We are able to place constraints on the properties of the molecular gas by studying the ratios of CO(1−0), (3−2), and (4−3). The molecular gas line ratios are defined as

$$\begin{eqnarray}&&{R}_{{ij}}=\displaystyle \frac{{I}_{\mathrm{CO}({\rm{i}}-({\rm{i}}-1))}}{{I}_{\mathrm{CO}({\rm{j}}-({\rm{j}}-1))}}\displaystyle \frac{{j}^{2}}{{i}^{2}}.\end{eqnarray} \tag{ 5 }$$

Line ratios R_ij are equivalent to line brightness temperature ratios, and can be used to determine the level of excitation of CO (Solomon & Vanden Bout 2005; Dannerbauer et al. 2009). R_ij = 1 implies the relative excitation of two lines is consistent with thermal excitation (e.g., Daddi et al. 2015). Based on the fluxes in Table 2, we find that for MACS 1931 as a whole, R₃₁ = 0.93 ± 0.16, and R₄₁ = 0.61 ± 0.10. R₃₁ in MACS 1931 is similar to the R₃₁ ∼ 0.8 values commonly observed in low-redshift BCGs (e.g., Edge 2001; Russell et al. 2016; Vantyghem et al. 2017).

The relative excitation levels of CO(3−2) and CO(4−3) in MACS 1931 places the excitation level of the molecular gas in this system similar to gas typically seen in low-redshift ULIRGs based on R₃₁, while R₄₁ is below the average for the local ULIRGs (Daddi et al. 2015). The CO spectral line energy distributions (SLEDs) for MACS 1931 and other astrophysical sources are shown in Figure 8. The MACS 1931 values are the red diamonds. The solid black line is the CO SLED for a star-forming galaxy predicted by Narayanan & Krumholz (2014), with CO excited by star formation with a surface density of ΣSFR = 1 M_⊙ yr⁻¹ kpc⁻², adapted from Daddi et al. (2015). For reference, the UV-derived ΣSFR of the MACS 1931 BCG is 0.83 ± 0.06 M_⊙ yr⁻¹ kpc⁻² (Fogarty et al. 2015). The submillimeter galaxy (SMG; blue points) and local ULIRG averages (gray points) are derived from the samples studied in Daddi et al. (2015) (see also Papadopoulos et al. 2012 and Bothwell et al. 2013), while the quasar-stellar object (QSO) average (purple triangles) is taken from Carilli & Walter (2013). For our estimate of the BCG CO(3−2) average flux (green point) we used sources in the literature with reported ALMA CO(1−0) and CO(3−2) observations, and adopted the scatter in these observations for the uncertainty. These include RXJ 1504-0248 (Vantyghem et al. 2018), 2A 0335 + 096 (Vantyghem et al. 2016), PKS0745-191 (Russell et al. 2016), and A1664 (Russell et al. 2014).

Zoom In Zoom Out Reset image size

We compare MACS 1931 to SMGs and ULIRGs since both have several characteristics in common with the MACS 1931 BCG. As with SMGs, the BCG is at the high end of the galactic stellar mass distribution and yet the observed large SFR contributes, at most, a few percent to their overall stellar mass budgets (Michałowski et al. 2010; González et al. 2011). SMGs are hypothesized to be the progenitors of elliptical galaxies and are often observed undergoing massive bursts of star formation as large as several × 1000 M_⊙ yr⁻¹ (Pope et al. 2006). Meanwhile, the IR luminosity of MACS 1931 is (1.4 ± 0.2) × 10¹² L_⊙, consistent with the definition of a ULIRG (Santos et al. 2016). While the majority of ULIRGs are fueled by galaxy interactions, the MACS 1931 BCG is typical of a ULIRG in terms of hosting a powerful dusty starburst and AGN (Lonsdale et al. 2006; Santos et al. 2016; Fogarty et al. 2017).

Figure 8 shows that there is mild tension between both the observed CO(3−2) and CO(4−3) transitions and the average behavior of the ULIRG sample. However, the MACS 1931 CO SLED is not consistent with the SMG sample. The SLED is also inconsistent with the Narayanan & Krumholz (2014) prediction for molecular gas excited by star formation with a similar ΣSFR to MACS 1931, although the CO(4−3) transition for MACS 1931 lies between the ULIRG SLED and the Narayanan & Krumholz (2014) prediction.

The tension between the MACS 1931 BCG and the ULIRG sample by itself does not imply a significant discrepancy between the BCG and the ULIRG population, since the CO SLED we plot is an average with considerable variation. However, it is noteworthy that, while R₃₁ implies that CO(3−2) is consistent with 1, implying extreme excitation conditions, R₄₁ is substantially less than this value. Papadopoulos et al. (2012) and Rosenberg et al. (2015) note that the conditions required for the levels of molecular gas excitation that produce R_j1 ≳ 1 at the scales of individual molecular clouds exist near extremely UV-bright young stars and supernova shocks, but producing these levels of excitation for an entire ULIRG requires either a very strong X-ray-dominated region from a QSO or large-scale mechanical turbulence or cosmic ray heating.

While AGN emission plays a role in the MACS 1931 CO SLED and provides a compelling explanation for the observed CO(3−2) line strength, QSO heating does not explain the average CO(4−3) level we observe. On the other hand, recent mechanical-mode AGN feedback can potentially provide both cosmic ray and turbulent heating to excite the MACS 1931 CO SLED. Comparison of the CO SLEDs suggest that the molecular gas excitation budget in MACS 1931 may be driven by a combination of star formation and processes that, where more strongly present, can drive ULIRG CO excitation ratios to >1 (Papadopoulos et al. 2012).

A related possibility is that some of the gas in MACS 1931 may be undergoing processes similar to what is observed in NGC 1275. Young stars are found to be associated with Hα and molecular filaments in NGC 1275, the central galaxy of Perseus, but there is often an offset of 0.5–1 kpc between young stars and filaments (Canning et al. 2014; Li et al. 2018). There are also many filaments in Perseus that are not actively forming stars (Salomé et al. 2008; Lim et al. 2017). These filaments may experience pressure support from turbulence and magnetic fields, and heating via infiltration and radiation from the surrounding ICM that suppresses the collapse of this gas into star-forming cores (Conselice et al. 2001; Fabian et al. 2008, 2011; Salomé et al. 2008). This process may help to drive CO(3−2) excitation but around the denser gas traced by CO(4−3) would potentially need to be weaker, and therefore provide an explanation for the CO SLED in MACS 1931.

In order to study the CO line ratios in MACS 1931 in more detail, in Figure 9 we map R₃₁ in all pixels where the value of ${R}_{31}/{\sigma }_{{R}_{31}}\geqslant 3$, after accounting for both rms uncertainties in the CO(1−0) and CO(3−2) images and the uncertainty in the absolute flux calibration. To examine R₃₁ in as much of the BCG as possible, we overlaid the ratio map obtained with beam-matched, naturally weighted images onto the ratio map obtained with beam-matched, 15 uv-tapered images. We superimposed the contours from the HST F336W UV image on the R₃₁ map.

Zoom In Zoom Out Reset image size

In the BCG core, there are several notable bright UV knots to the north of the AGN and a depression in UV output to the south, which we attribute to a dust lane. R₃₁ loosely traces both sets of features, suggesting local enhancement of the CO lines ratios by star formation to the north and by the process involved in creating the dust lane to the south. R₃₁ is lower in the tail than in the core, and roughly drops off with increasing distance from the AGN point source. However, the tail and regions to the south and east of the core show areas of enhancement in R₃₁.

Spatial variation in the excitation of CO in MACS 1931 has some similarities to the variation observed by Lim et al. (2017) in NGC 1275, who found that R₃₂ varied between ∼1 near the center of the galaxy and ∼0.5 further out. Furthermore, R₂₁ values observed in NGC 1275 by Salomé et al. (2008) varied between 0.6 and 1.7 in the center, but were observed to be as high as 2.4 at offset positions. MACS 1931 similarly shows enhancements of R₃₁ at positions several 10 s of kpc from the BCG center.

There are likely varying contributions to the molecular gas excitation in cool-core BCGs from star formation, from the AGN, and from collisional excitation due to AGN feedback moderated turbulence (e.g., Lim et al. 2017). This variation dominates the error bar in the average R₃₁ value for the BCG subsample (green data point) shown in Figure 8.

Molecular gas excitation is also a potentially useful way to compare BCGs in cool-core clusters like MACS 1931 to BCG progenitors in high-redshift protoclusters. Protoclusters can contain galaxies with even more massive molecular gas reservoirs than that found in MACS 1931, but the origins of the molecular gas and how it interacts with the cluster environment may be substantially different (Emonts et al. 2016; Dannerbauer et al. 2017). For example, compared to MACS 1931, the molecular gas feature in the central radio galaxy of a z = 2.2 protocluster known as the Spiderweb Galaxy is significantly more excited, and contains 10 times as much molecular gas ((2 ± 0.2) × 10¹¹ M_⊙) (Emonts et al. 2018). R₄₁ in the Spiderweb is ∼1, suggesting more extreme molecular gas excitation conditions than in MACS 1931. However, among the studied high-redshift radio galaxies (a population which is also thought to contain progenitors of low-redshift BCGs), a wide range of molecular gas masses (∼10¹⁰–10¹¹ M_⊙) have been reported, along with a wide range of excitation conditions (Miley & De Breuck 2008). The significant differences in both the amount of gas and excitation of the gas suggest that the processes involved in producing star formation and molecular gas may differ in cool-core clusters with quasar-like AGNs and in protoclusters.

Observations of dust continuum emission in Bands 6 and 7 reveal that BCG dust is concentrated in regions that also contain multiphase gas and star formation. The dust-to-molecular gas ratio in the MACS 1931 BCG is ∼0.01–0.1, based on the total molecular gas mass reported here and the dust mass derived in Fogarty et al. (2017). A similar dust-to-gas ratio in the range ∼0.04–0.1 is observed in the Phoenix cluster BCG (Russell et al. 2017a; Fogarty et al. 2017). However, the relative distributions of dust and gas in MACS 1931 do not appear to be uniform, since the relative far-IR emission in Bands 6 and 7 due to dust along the Hα tail is much greater than that in seen in the core (excluding possible contributions to dust emission from the AGN point source that were subtracted out).

There is substantial variation in the relative far-IR flux due to dust compared to CO, UV, and Hα flux in the Hα tail extending northwest of the BCG. In Figure 10, we plot the normalized far-IR, CO(3−2), UV, and Hα+[N ii] surface brightness profiles in a 3'' wide rectangular slice along the Hα tail. Dust emission is elevated in the ∼10 kpc of the tail furthest from the core compared to the UV, Hα+[N ii], and CO(3−2) surface brightness profiles. This result suggests that the dust concentration is enhanced in the outskirts of the tail.

Zoom In Zoom Out Reset image size

Dust in the BCG either forms in situ, is launched out of the BCG core by jets, or is the result of a combination of processes. In situ cooling requires that dusty winds from AGB stars be preserved long enough to end up in the surrounding dusty nebular structure (Li et al. 2019). This source of dust may be supplemented by supernova production boosted by the starburst (Dunne et al. 2003). Alternatively, dust in the BCG can arise from denser concentrations of dusty material dredged out of the center of the BCG by AGN jets (Fabian et al. 1994). In this scenario, dusty material seeds the low-entropy gas uplifted by jets, which then induces condensation of the dust out of the resulting molecular clouds (Fabian et al. 1994; Voit & Donahue 2011; Voit et al. 2015). The observed distribution of the dust in MACS 1931 relative to the atomic and molecular gas and UV emission from young stars suggests some uplift-dominated process is taking place, with dust emission most prominent in the core and at the furthest extremity of the tail, which is presumably the most recently uplifted body of material. In all likelihood, these processes co-mingle, with seed grains supplied by both jets and the AGB population stimulating further condensation of dust out of the molecular gas.

Like the origin of dust, the persistence of dust in cool-core BCGs is a topic of debate, since sputtering induced by interactions with the hot ICM is expected to destroy dust on timescales of ≲10⁶ yr (Voit & Donahue 2011). However, the very cold, ∼10 K, temperatures that we derive for some of the dust in MACS 1931 suggests that a fraction of the dust in the BCG must be protected from sputtering or is in pockets of the ICM that are substantially colder than the 4.78 ± 0.64 keV temperature of the ICM plasma around the multiphase tail observed in Ehlert et al. (2011).

Montier & Giard (2004) studied collisional heating of dust by the ICM, and concluded that, for a wide range of ICM densities and temperatures above ∼1.7 keV, the equilibrium temperature for dust interacting with the ICM does not depend substantially on grain size and is a function ICM density, n_e. Assuming that n_e > 0.1 cm⁻³ (corresponding to the radial profile bin between ∼10 and ∼30 kpc in Ehlert et al. 2011), the relation Montier & Giard (2004) derive implies the dust equilibrium temperature must be at least 38 K. A similar analysis by Dwek et al. (1990) provides the same conclusion.

Meanwhile, Popescu et al. (2000) found that for small grain sizes the dust temperature distribution may peak at ∼10 K, but they studied the core of the Virgo Cluster, which has an ICM temperature of ∼1.1 keV. In MACS 1931, a plausible scenario may thus be that the ∼10 K dust in the ICM is either interacting with local volumes of ≲1 keV material and has a large enough proportion of small grains that this population determines the Band 6 and 7 fluxes, or this dust is somehow protected from interacting with the hot ICM.

In the former case, we may hypothesize that the cold dust itself is a potential coolant. While the role of dust in promoting ICM cooling is still a subject of debate, recent findings suggest dust is capable of enhancing ICM cooling by up to 50% (Vogelsberger et al. 2018). The multiphase tail to the northwest of the MACS 1931 BCG is located on a prominence in both ICM metallicity and X-ray flux, and the ICM surrounding the tail has the lowest entropy (13–16 keV cm²) of anywhere in the cluster (Ehlert et al. 2011). The extremely low entropy in the ICM immediately surrounding the tail (see Figure 5(d) in Ehlert et al. 2011) suggests the potential for further condensation, although it may also imply a conductive or turbulent mixing layer. In the event that the very cold dust in the multiphase tail is tracing pockets of ICM that are several keV colder than this already low-entropy gas, it is possible that the dust itself is promoting condensation of the ICM onto the multiphase tail, and the growth of dust out of the molecular gas is sufficient to allow the extended dust to persist.

As previously observed in Giacintucci et al. (2014) and Yu et al. (2018), there is a radio mini-halo in the core of MACS 1931 which may be associated with the feedback processes occurring in this system. Radio mini-halos are small, typically diffuse, structures around cool-core BCGs that extend across some or all of the cluster core (Gitti 2016). While radio mini-halos are thought to occur in most cool-core cluster cores (Giacintucci et al. 2017), the mini-halo in MACS 1931 is unusual for having two compact 1.5 GHz sources (Giacintucci et al. 2014): an 11.6 mJy source on the position of the BCG core and a 2.5 mJy source on the position of the tail as shown Figure 11. We note that the 1.5 GHz flux predicted by the Schmitt et al. (2006) relation for an SFR of 250 M_⊙ yr⁻¹ is ∼1 mJy, so neither component can be easily explained solely by star formation in the BCG.

Zoom In Zoom Out Reset image size

While the origins of radio mini-halos in cool-core clusters are not well known, observations and simulations suggest that these objects are due to turbulent re-acceleration of relativistic electrons in the ICM, which arise from the relativistic plasma injected into the system by AGN jets (Fujita et al. 2007; ZuHone et al. 2013; van Weeren et al. 2019). The presence of the smaller compact 1.5 GHz component is striking because of its correspondence with the multiphase tail. It is likely that this feature arises out of interactions between the relativistic electron population of the ICM and one or more components of the multiphase tail, although it is not clear what is driving the acceleration of the electron population. We point out this correlation because it is intriguing, and merits further investigation.