Testing the Limits of AGN Feedback and the Onset of Thermal Instability in the Most Rapidly Star-forming Brightest Cluster Galaxies

Our "extreme cooling sample" originally consisted of six of the most star-forming (SFR ≳ 100 M_⊙yr⁻¹) and strongly cooling (${\dot{M}}_{\mathrm{cool}}\gtrsim 1000$ M_⊙yr⁻¹) BCGs studied in McDonald et al. (2018). These systems include the Phoenix cluster (z = 0.596), H1821+643 (z = 0.299, hereafter "H1821"), IRAS 09104+4109 (z = 0.442, "IRAS09104"), Abell 1835 (z = 0.252, "A1835"), RX J1532.9+3021 (z = 0.363, "RXJ1532"), and MACS 1931.8-2634 (z = 0.352, "MACS1931"). SFRs for these systems, presented in McDonald et al. (2018), were aggregated from multiple literature sources that utilize various SFR measures (e.g., Hα, UV, or IR flux). With the exception of RXJ1532 and MACS1931, with 16+ band HST imaging from CLASH ¹³ allowing careful stellar population modeling (Fogarty et al. 2017), AGN contamination was thought to be affecting the SFR measurements of most of these systems, making the quoted SFRs less secure, even after spatial or spectral decomposition attempts. The new observations we present below will provide more secure SFRs for these systems. In addition, RBS797 (z = 0.354), which was previously estimated to have a low SFR relative to its ICM cooling rate (${\dot{M}}_{\mathrm{cool}}\,\sim 1000$ M_⊙yr⁻¹), was also added to our sample upon measuring a much higher updated SFR for it in a new, separate HST narrowband program, as we will show below in Section 3.2.

We obtained new HST optical data for our sample during programs GO15661 and GO16494 (PI: McDonald), and GO16001 (PI: Gitti). These included narrowband data using the ramp filters on the Advanced Camera for Surveys (ACS)/Wide Field Camera (WFC). These data were taken over 40 orbits, to which we include an additional 13 orbits on the Phoenix cluster (GO15315) that were previously published (McDonald et al. 2019), as well as two orbits each from the CLASH survey (Postman et al. 2012) for MACS1931 (GO12456) and RXJ1532 (GO12454), and less than one additional orbit for RBS797 from SNAP program 10875 (PI: Ebeling). Our observing strategy was designed to tune the ramp filters to be centered on the [OII]λ λ3726, 3729 doublet at the redshift of each cluster. This choice also avoids contamination from any other strong cooling lines to facilitate interpretation as gas that will eventually form stars. The broadband filters were chosen in such a way as to be free of contamination from the strongest expected emission lines ([OII], [OIII]λ λ4959, 5007, and Hα) to better model and subtract the underlying continuum from the narrowband observations (described in more detail below in Section 2.2.1). A summary of the new and archival data used in this study can be found in Table 1, where we list the exposure times, proposal IDs, and filters used for each source, including the wavelength each ramp filter was tuned to for observing redshifted [OII] emission. All of the HST data used in this paper can be found in MAST:10.17909/vdq6-9693.

For each system, we used STScI's DrizzlePac software package ¹⁴ to further process the individual calibrated, flat-fielded exposures. For all visits of a single filter, we used SExtractor (Bertin & Arnouts 1996) to create point-source catalogs that were passed into TweakReg (v1.4.7) in order to register all exposures to the same World Coordinate System at the subpixel level. We then used AstroDrizzle (v3.1.6) to remove geometric distortions, correct for sky background variations, flag cosmic rays, and combine individual exposures. This procedure was repeated for both of the broadband filters used for each target. The ramp filter exposures of each system were combined in a similar manner, except with an additional initial step of running AstroDrizzle first with only cosmic-ray flagging in order to be able to run SExtractor, as the raw exposures are usually dominated by cosmic rays due to the filters' low throughputs and narrow fields of view (see ACS instrument science reports ¹⁵ ; Lucas & Hilbert 2015; Lucas 2015). After registering, cleaning and combining the exposures for each broadband and ramp filter separately, we create point-source catalogs and run TweakReg again on each of the now combined filter data to align the images across all filters.

2.2.1. Continuum Subtraction

In order to do the continuum subtraction, we first compose a spectral energy distribution (SED) for each pixel. This is done by convolving the broadband filter bandpasses ${{ \mathcal T }}_{{\rm{X}}}$ (obtained with pysynphot) used for each cluster with both a redshifted young (10 Myr) and old (6 Gyr) stellar population spectral template obtained with STARBURST99 (Leitherer et al. 1999). The predicted flux from these convolutions is then scaled to match the observed flux in each broadband observation. Because the number of convolved stellar templates is equal to the number of broadband filters, a unique linear combination exists that allows us to convert between blue and red bandpass fluxes (F_B and F_R) to "young" (${ \mathcal Y }$) and "old" (${ \mathcal O }$) template fluxes:

$$\begin{eqnarray}&&\left[\begin{array}{l}{F}_{{\rm{B}}}\\ {F}_{{\rm{R}}}\end{array}\right]=\left[\begin{array}{ll}\displaystyle \frac{\int { \mathcal Y }\ast {{ \mathcal T }}_{{\rm{B}}}d\lambda }{\int {{ \mathcal T }}_{{\rm{B}}}d\lambda } & \displaystyle \frac{\int { \mathcal O }\ast {{ \mathcal T }}_{{\rm{B}}}d\lambda }{\int {{ \mathcal T }}_{{\rm{B}}}d\lambda }\\ \displaystyle \frac{\int { \mathcal Y }\ast {{ \mathcal T }}_{{\rm{R}}}d\lambda }{\int {{ \mathcal T }}_{{\rm{R}}}d\lambda } & \displaystyle \frac{\int { \mathcal O }\ast {{ \mathcal T }}_{{\rm{R}}}d\lambda }{\int {{ \mathcal T }}_{{\rm{R}}}d\lambda }\end{array}\right]\left[\begin{array}{l}{c}_{{ \mathcal Y }}\\ {c}_{{ \mathcal O }}\end{array}\right].\end{eqnarray} \tag{ 1 }$$

This equation can be inverted to solve for the vector ${\left[\begin{array}{ll}{c}_{{ \mathcal Y }} & {c}_{{ \mathcal O }}\end{array}\right]}^{T}$, the coefficients by which to scale the template spectra to match the observed blue and red broadband fluxes. Because the templates are defined over a broad range of wavelengths, we can integrate the scaled composite young plus old templates over the ramp filter bandpass to estimate the expected narrowband continuum flux (F_N ):

$$\begin{eqnarray}\begin{array}{rcl}{F}_{N} & = & \left[\begin{array}{ll}\displaystyle \frac{\int { \mathcal Y }\ast {{ \mathcal T }}_{N}d\lambda }{\int {{ \mathcal T }}_{N}d\lambda } & \displaystyle \frac{\int { \mathcal O }\ast {{ \mathcal T }}_{N}d\lambda }{\int {{ \mathcal T }}_{N}d\lambda }\end{array}\right]\left[\begin{array}{l}{c}_{{ \mathcal Y }}\\ {c}_{{ \mathcal O }}\end{array}\right]\\ & = & \displaystyle \frac{\int \left({c}_{{ \mathcal Y }}\cdot { \mathcal Y }+{c}_{{ \mathcal O }}\cdot { \mathcal O }\right)\ast {{ \mathcal T }}_{N}d\lambda }{\int {{ \mathcal T }}_{N}d\lambda }.\end{array}\end{eqnarray} \tag{ 2 }$$

Finally, we can subtract this predicted narrowband flux from the observed narrowband flux to get a continuum-subtracted [OII]-only flux. This SED fitting and continuum subtraction is done pixel-by-pixel to get [OII]-only maps for each of the clusters in our sample. The procedure is illustrated in Figure 1.

Zoom In Zoom Out Reset image size

Archival Chandra X-ray data were used for each of the clusters in our sample in order to qualitatively compare to the features seen in the HST [OII] maps. A summary of the observations used for this analysis is presented in Table 2. All Chandra observations were reduced using the Chandra Interactive Analysis of Observations (CIAO) v4.10 software package, along with CALDB v4.8.0. The observations were reprocessed using chandra_repro, applying the latest gain and charge transfer inefficiency corrections, with improved background filtering applied to those observations taken in the VFAINT telemetry mode. Observations from multiple ObsIDs were reprojected, exposure-corrected, and combined using the merge_obs routine, over an energy range of 0.5–7 keV. Point sources were removed via a wavelet decomposition of the merged images using the wavdetect script, after which periods of high background were excluded using lc_clean. Blank-sky background files used for background subtraction were obtained using the CIAO blanksky script. These blank-sky files were renormalized to have the same high-energy particle rate as the observations over the 9.5–12 keV range, where Chandra's effective area is low enough so that any flux is mostly due to the particle background, using the blanksky_sample script.

2.3.1. Thermodynamic Profiles

For our spectral analysis, spectra were extracted from concentric annuli centered on the X-ray peak in each system, using broad and fine bin widths, as in, e.g., McDonald et al. (2019). The broader annuli were chosen to contain ≳2000 counts to permit precise temperature measurements. The spectra were extracted from the observation and blank-sky background files over the energy range of 0.5–7 keV using specextract, with an identical off-source region also being extracted for both. The spectra from multiple ObsIDs in each bin were fit with the XSPEC v12.10.1 spectral fitting package (Arnaud 1996) using the Cash statistic (Cash 1979) and Levenberg–Marquardt algorithm. Each on-source spectrum was fit using a combined spectral model PHABS∗(APEC + BREMSS) + APEC, where the first APEC component is to model thermal emission from the optically thin plasma in the ICM (Smith et al. 2001) along with Galactic photoelectric absorption using PHABS. We further model the background with a second APEC component (fixed at kT = 0.18 keV, Z = Z_⊙, z = 0) to model soft Galactic X-ray emission, as well as a BREMSS model (fixed at kT = 40 keV) to model unresolved point sources (e.g., McDonald et al. 2019). These two background components are jointly fit with the off-source spectra, which are only modeled with PHABS*BREMSS + APEC, with normalizations scaled by the extraction area of the on-source spectra. Galactic Hi column densities (N_H) for each system were taken from Kalberla et al. (2005), redshifts were fixed to their literature values, and metallicities for the on-source APEC components were fixed at Z = 0.3Z_⊙. Only APEC temperatures and normalizations were allowed to vary.

The resulting coarse temperature profiles extracted along each annular bin were fit with the analytical model of Vikhlinin et al. (2006):

$$\begin{eqnarray}&&{T}_{3{\rm{D}}}(r)={T}_{0}\displaystyle \frac{{\left(r/{r}_{\mathrm{cool}}\right)}^{\alpha }+{T}_{\min }/{T}_{0}}{{\left(r/{r}_{\mathrm{cool}}\right)}^{\alpha }+1}\displaystyle \frac{{\left(r/{r}_{t}\right)}^{a}}{{[1+{\left(r/{r}_{t}\right)}^{b}]}^{c/b}},\end{eqnarray} \tag{ 3 }$$

where a, b, and c model the outer regions of the profile with a flexible broken power law. The profile transitions to an inner cool-core component at around r_t , with the cool core defined by T₀, ${T}_{\min }$, r_cool, and α. All fitting parameters are initialized to the average values found in Vikhlinin et al. (2006). This 3D temperature was projected along the line of sight and over the width of each annular bin, and then fit to the observed profile using the python package LMFIT (Newville et al. 2014). The corresponding best-fit temperature profile was then interpolated along a finer radial binning from which we extracted another profile. In these finer spectral fits, we keep the temperatures fixed at these interpolated values. Only the APEC normalization is allowed to vary in order to reduce the degrees of freedom and more precisely constrain the density profile. The APEC normalization η was converted to an emission measure (EM) profile using EM(r) $\equiv \int {n}_{{\rm{e}}}{n}_{{\rm{p}}}{dV}=\eta \times 4\pi \times {10}^{14}{\left[{D}_{{\rm{A}}}(1+z)\right]}^{2}$, where n_e and n_p are the electron and proton number densities, and D_A is the angular diameter distance at the cluster redshift. The EM profile was fit with the following model from Vikhlinin et al. (2006):

$$\begin{eqnarray}\begin{array}{rcl}{n}_{{\rm{e}}}{n}_{{\rm{p}}}(r) & = & {n}_{0}^{2}\displaystyle \frac{{\left(r/{r}_{{\rm{c}}}\right)}^{-\alpha }}{{\left[1+{\left(r/{r}_{{\rm{c}}}\right)}^{2}\right]}^{3\beta -\alpha /2}}\displaystyle \frac{1}{{\left[1+{\left(r/{r}_{{\rm{s}}}\right)}^{\gamma }\right]}^{\epsilon /\gamma }}\\ & & +\displaystyle \frac{{n}_{02}^{2}}{{\left[1+{\left(r/{r}_{{\rm{c}}2}\right)}^{2}\right]}^{3{\beta }_{2}}},\end{array}\end{eqnarray} \tag{ 4 }$$

which is a modified double-beta model with a cusp, rather than a flat core (defined by n₀, r_c, α, β), a steeper outer profile slope (defined by r_s, γ, and ), and a cool-core component (defined by n₀₂, r_c2, and β₂). In our fits, γ = 3 remains fixed, and all other parameters are allowed to vary and initialized to typical parameters found in Vikhlinin et al. (2006). This 3D model was also projected along the line of sight. From this fit, the electron number density was calculated assuming abundances from Anders & Grevesse (1989) for a fully ionized 0.3Z_⊙ abundance plasma, such that n_e/n_p = 1.2.

Additional thermodynamic profiles are derived analytically from the fitted density and temperature profiles. We calculate profiles for the total pressure (P = (n_e + n_p)k_B T), pseudo-entropy ($K={k}_{{\rm{B}}}{{Tn}}_{{\rm{e}}}^{-2/3}$), cooling time (${t}_{\mathrm{cool}}=\tfrac{3}{2}\tfrac{({n}_{{\rm{e}}}+{n}_{{\rm{p}}}){k}_{{\rm{B}}}T}{{n}_{{\rm{e}}}{n}_{{\rm{p}}}{\rm{\Lambda }}({k}_{{\rm{B}}}T,Z)}$), and freefall time (${t}_{\mathrm{ff}}=\sqrt{2r/g}$). We used the cooling function Λ(k_B T, Z) from Sutherland & Dopita (1993), as parameterized by Tozzi & Norman (2001; see also Guo & Oh 2008; Parrish et al. 2009). The gravitational acceleration used in calculating the freefall time was obtained by modeling the cluster potential with the sum of an isothermal sphere at small radii and a Navarro–Frenk–White profile (Navarro et al. 1997) at larger radii, assuming a velocity dispersion of σ_v = 350 km s⁻¹ typical of BCGs (e.g., Von Der Linden et al. 2007; Sohn et al. 2020). We compare the profiles of each of our sources to measurements from the literature to confirm agreement between our thermodynamic modeling. Literature sources we referenced for comparison include: McDonald et al. (2019) for Phoenix, Russell et al. (2010) for H1821, O'Sullivan et al. (2012) for IRAS09104, Cavagnolo et al. (2009) for A1835, RXJ1532 and MACS1931, as well as Ehlert et al. (2011) for MACS1931, and Cavagnolo et al. (2011) and Doria et al. (2012) for RBS797. Data for the Phoenix cluster and H1821+643 were obtained from M. McDonald and H. Russell, respectively (via private communication), and profiles were not extracted independently here, only modeled analytically, due to the complicated nature of the AGN's contribution in these systems.

Our new, reduced broad- and narrowband HST observations can be found in Figure 2. The continuum-subtracted [OII] maps in the middle panels reveal with remarkable detail the intricate morphology of the warm (T ∼ 10⁴ K) ionized gas within each of the BCGs in our sample that will eventually form stars. These maps represent the highest angular resolution view of the star-forming material at optical wavelengths in each of these strongly cooling clusters to date (with the exception of Phoenix, which was already presented in McDonald et al. 2019). Each of these nebulae are extended on the order of tens of kiloparsecs in projection, ranging from R_[OII] = 23–60 kpc, making them among the most extended emission-line nebulae (see, e.g., McDonald et al. 2010; Hamer et al. 2016). The most extended of these filaments is found in Phoenix, where the [OII] filaments reach up to 60 kpc (McDonald et al. 2019), similar to the well-studied Hα filament system found in the extremely deep observations of NGC1275 at the center of the nearby (z = 0.018) Perseus cluster (Conselice et al. 2001; Fabian et al. 2003), with R_Hα ≈ 40 kpc.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The morphologies of the [OII] maps appear to be extending either in various directions (as is the case in Phoenix, A1835, and RXJ1532) or in predominantly one direction (as in H1821, IRAS09104, MACS1931, and RBS797). The cause of these different morphologies may be investigated by looking to observations from other wavelengths. For instance, in the Phoenix cluster, the X-shaped morphology of the pairs of filaments to the north and south are coincident with the rims of X-ray cavities carved out by radio bubbles, as detailed in McDonald et al. (2019; see Figure 2). Cold molecular gas coincident with Hα filaments, as measured by ALMA (Russell et al. 2017a), suggest that in a few systems buoyantly rising radio bubbles are promoting cooling in their wake as they climb out of their deep cluster potentials in a few systems (e.g., Revaz et al. 2008; Pope et al. 2010; McNamara et al. 2016). This picture has been similarly suggested as the mechanism for cooling in A1835, where the [OII] filaments presented here for the first time are also coincident with molecular gas rising behind the location of X-ray cavities (McNamara et al. 2014). For each of the clusters in our sample, we can see the morphology of the [OII] overlaid on top of the X-ray maps of the cluster cores, along with the location of known X-ray cavities from the literature in Figure 2. Cool gas in some clusters appears to be extended in mostly one direction, in some cases completely perpendicular to the direction of the known X-ray cavities, as in MACS1931, where there is also molecular gas that traces the [OII] and Hα + [NII] morphology (Fogarty et al. 2019). A possible explanation for the triggering of star formation that is not stimulated by rising radio bubbles could be a recent gas-rich merger that deposited the cooling material or created turbulence to promote existing material to cool as per the CCA model. Clearly, the location of bubbles does not necessarily predict the direction of cooling. Furthermore, even in the case of Phoenix, where there is obvious agreement between the position angles of cavities and star-forming filaments, the [OII] filaments extend even beyond the maximum radius of the cavities, implying that radio bubbles do not tell the whole story. We will discuss the implications of the presence of bubbles on the development of thermal instabilities further in Section 4.2.1.

One of the main goals of our new HST observations was to secure more precise SFRs by being able to more faithfully isolate the morphology of star-forming filaments and exclude likely regions of high AGN contamination. To that end, we extract an [OII] flux f_[OII] from polygonal apertures designed to encompass all of the flux while maximizing the signal-to-noise for each of the [OII] nebulae. We applied an intrinsic extinction correction to each of our rest-frame [OII] flux measurements based on the E(B − V) measurements by Crawford et al. (1999), assuming R_V = 3.1. In the case of Phoenix, we used an E(B − V) =0.24 ± 0.10 from McDonald et al. (2014). Crawford et al. (1999) measured specific reddening values for A1835 of E(B − V) = 0.38 ± 0.04 and for RXJ1532 of E(B − V) =0.21 ± 0.03. For the rest of our sample, we take the distribution of reddening values from the entire sample of Crawford et al. (1999; Table 4, Column 7), and use the first, second, and third quartiles to apply a reddening of $E(B-V)={0.27}_{-0.07}^{+0.16}$. Following this intrinsic extinction correction, we also apply a Galactic extinction correction at the redshifted [OII] wavelengths using the values listed in the NASA/IPAC Extragalactic Database, which are based on Sloan $r^{\prime}$ band data from Schlafly & Finkbeiner (2011). Following these corrections, we converted our corrected [OII] fluxes to SFRs using the scaling relations found in Kewley et al. (2004):

$$\begin{eqnarray}&&{\mathrm{SFR}}_{[{\rm{O}}{\rm\small{II}}]}=6.58\times {10}^{-42}{L}_{[{\rm{O}}{\rm\small{II}}]}\,\,{M}_{\odot }\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 5 }$$

where L_[OII] = f_[OII] · 4π D_L(z)² is in units of erg s⁻¹, and D_L (z) is the redshift-dependent luminosity distance. These new [OII] SFRs can be found in Table 3. We find good agreement between our new measurements and the literature, particularly in the case of Phoenix (McDonald et al. 2019), as well as MACS1931 and RXJ1532 (Fogarty et al. 2017), IRAS09104 (Ruiz et al. 2013), and A1835. Previous SFR values for RBS797 based on UV fluxes (Cavagnolo et al. 2011) are much lower than our [OII] SFR measurements, likely due to no extinction correction being made.

Table 3. New [OII] SFRs for Our Sample

Name	${\mathrm{log}}_{10}{\mathrm{SFR}}_{[{\rm{O}}{\rm\small{II}}]}$	${\mathrm{log}}_{10}{\dot{M}}_{\mathrm{cool}}$
	(M⊙ yr−1)	(M⊙ yr−1)
Phoenix	2.76 ± 0.20	3.23 ± 0.08
H1821+643	2.08 ± 0.30	2.65 ± 0.05
IRAS 09104+4109	2.30 ± 0.30	3.01 ± 0.04
Abell 1835	2.05 ± 0.08	3.07 ± 0.06
RX J1532.9+3021	1.95 ± 0.06	3.03 ± 0.04
MACS 1931.8-2634	2.19 ± 0.30	3.03 ± 0.10
RBS 797	1.69 ± 0.30	3.15 ± 0.24

Download table as:  ASCIITypeset image

Special care had to be taken to extract an accurate f_[OII] for H1821, as the bright quasar produced strong negative residuals in the continuum subtraction due to the diffraction spikes, as seen in Figure 2. These negative residuals were masked in our extraction region. Also, we might expect that some of the ionized [OII] emission is due to radiation from the quasar and not from star formation. To account for this, we initially tried removing a central region up to a radius defined by the Strömgren sphere:

$$\begin{eqnarray}&&{r}_{\mathrm{Str}\ddot{{\rm{o}}}\mathrm{mgren}}={\left[\displaystyle \frac{3{N}_{{\rm{i}}}}{4\pi {\alpha }_{{\rm{B}}}\langle {n}_{{\rm{e}}}{\rangle }^{2}}\right]}^{1/3},\end{eqnarray} \tag{ 6 }$$

where α_B = 2.6 × 10⁻¹³ cm³ s⁻¹ is the coefficient for case B recombination (i.e., only net recombinations, thus excluding transitions directly to the ground state), 〈n_e〉 = 300 cm⁻³ is the average electron density of cool line-emitting gas in cluster cores (e.g., McDonald et al. 2012a), and N_i is the ionizing photon rate. The ionizing photon rate was estimated by finding a best-fit relation of $\mathrm{log}{N}_{{\rm{i}}}=1.05\mathrm{log}{L}_{\mathrm{bol}}+7.46$ between the bolometric luminosity L_bol and ionizing photon rate for the quasar sample in Elvis et al. (1994; see their Table 14). Using a bolometric luminosity of L_bol = 2 × 10⁴⁷ erg s⁻¹ for H1821 (Russell et al. 2010), we find an approximate ionizing photon rate of N_i ≈ 10⁵⁷ s⁻¹ and a corresponding Strömgren sphere radius of about 0.7 kpc = 016. The amount of quasar contamination to the [OII] flux within such a small radius is negligible, and even more so for the rest of the systems in our sample. However, given that this Strömgren sphere calculation presumes a homogeneous medium, and that the filling factor of the line-emitting gas is likely low, the effective Strömgren sphere radius should be larger. In addition, the Airy diffraction pattern induced by the bright quasar in H1821 also leaves behind positive residuals that cannot be included in the SFR calculation. Instead, in lieu of the complexity of modeling the HST PSF, we calculate a radius of 18 to exclude from the center of this source (red dashed circle in Figure 2), which corresponds to an encircled energy fraction of ∼99.8% of an on-axis point source. As a result, the SFR we measure for H1821 can be considered a lower limit, though the extent of the missing flux from young stars is uncertain, as much of the ionization from this masked region may come from the quasar.

In addition to our sample of seven clusters, we compare their SFRs to those of other cool-core clusters from the literature. In particular, we searched for systems with available Hα flux measurements of strongly cooling groups and clusters of galaxies, as Hα is a more readily available measurement for lower-redshift clusters in the current literature. We compiled the Hα fluxes/luminosities from Hamer et al. (2016), McDonald et al. (2010), Gomes et al. (2016), Buttiglione et al. (2009), Lakhchaura et al. (2018), Cavagnolo et al. (2009), Donahue et al. (1992), Crawford et al. (1999), Wang et al. (2010), and Heckman et al. (1989), prioritizing first the sources with tunable filter or integral-field spectroscopy, and then those with long-slit spectroscopy. For sources that quoted an Hα + [NII] flux, we converted to a pure Hα flux by using the ratio of Hα/(Hα + [NII]) = 0.76 (e.g., McDonald et al. 2010). We homogenized each of these literature flux measurements to conform to our same chosen cosmological parameters. In almost every case, a Galactic extinction correction had already been applied, after which we applied an intrinsic correction according to the distribution of E(B − V) values from Crawford et al. (1999), which becomes the largest source of uncertainty in each of these measurements. Finally, we attempt to remove the contamination to the Hα flux from evolved stars (e.g., planetary nebulae, supernova remnants, and AGB and HB stars). This contamination is likely to significantly contribute to the SFR especially for the weakly cooling clusters. To account for this contamination, we follow the same procedure described in McDonald et al. (2021), where more details can be found.

The SFRs for our extreme cooling sample, as well as the reference samples listed above, are plotted against the X-ray cooling rates (${\dot{M}}_{\mathrm{cool}};$ compiled by McDonald et al. 2018) in Figure 3, and can also be found in Tables 3 and 4, including the intermediate homogenization of the literature luminosity values. The classically inferred ICM cooling rates here are defined as ${\dot{M}}_{\mathrm{cool}}=\tfrac{{M}_{\mathrm{gas}}(r\lt {r}_{\mathrm{cool}})}{{t}_{\mathrm{cool}}}$, where r_cool is the radius where the cooling time is < 3 Gyr, and t_cool is 3 Gyr. This value of the cooling radius/timescale is chosen to more closely probe the cluster core where cooling actually occurs. There is typically good agreement between these estimates of cooling rates and luminosity-based rates, while other choices of the cooling radius (e.g., r_cool at t_cool = 7.7 Gyr) are essentially arbitrary but useful conventions for comparing to the literature. Spectroscopic-based cooling rates are more accurate measurements of how much gas is actually cooling, which is typically less than that inferred by "classical" (i.e., maximal) cooling rates and would tend to bring ${\dot{M}}_{\mathrm{cool}}$ measurements within an order of magnitude or less of the SFR measurements. However, spectroscopic cooling rates are much harder to measure and usually result in upper limits for cooling clusters. Our choice of cooling rate method is useful as an upper limit to the total rest mass of gas that could potentially cool, allowing energy budget considerations that can more conveniently reveal the impact of heating from AGN feedback.

Zoom In Zoom Out Reset image size

Table 4. Hα SFRs and ${\dot{M}}_{\mathrm{cool}}$ Data Used in Figure 3

Name	z	Hαlit	LHα,homog.	LHα,extinc.	LHα,corr.	${\mathrm{log}}_{10}{\mathrm{SFR}}_{{\rm{H}}\alpha }$	${\mathrm{log}}_{10}{\dot{M}}_{\mathrm{cool}}$
			1040erg s−1	1040erg s−1	1040erg s−1	log10 M⊙ yr−1	log10 M⊙ yr−1
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
2A 0335+096	0.035	8.3a	8.53	${15}_{-2}^{+7}$	${15}_{-2}^{+7}$	${0.16}_{-0.06}^{+0.16}$	2.26 ± 0.01
3C295	0.464	2.2e-15h	134.32	${240}_{-30}^{+100}$	${240}_{-30}^{+100}$	${1.4}_{-0.1}^{+0.2}$	2.98 ± 0.04
A0085	0.056	1.6a	1.64	${2.9}_{-0.4}^{+1.3}$	${2.3}_{-0.4}^{+1.3}$	−0.66${}_{-0.08}^{+0.19}$	1.94 ± 0.01
A0133	0.056	1.2a	1.23	${2.2}_{-0.3}^{+1.0}$	${1.8}_{-0.3}^{+1.0}$	−0.77${}_{-0.08}^{+0.19}$	1.79 ± 0.01
A0478	0.088	23a	23.57	${42}_{-6}^{+18}$	${42}_{-6}^{+18}$	${0.60}_{-0.06}^{+0.16}$	2.64 ± 0.01
A0496	0.033	3.1a	3.18	${5.7}_{-0.8}^{+2.5}$	${5.3}_{-0.8}^{+2.5}$	−0.29${}_{-0.07}^{+0.16}$	1.75 ± 0.03
A1204	0.171	59a	60.24	${110}_{-10}^{+50}$	${110}_{-10}^{+50}$	${1.0}_{-0.1}^{+0.2}$	2.6 ± 0.03
A1361	0.117	13.5d	6.89	${12}_{-2}^{+5}$	${12}_{-2}^{+5}$	${0.067}_{-0.062}^{+0.158}$	1.66 ± 0.19
A1413	0.143	−2.38j	−2.38	<2.38	<2.38	<−0.42	1.74 ± 0.12
A1650	0.084	0.03b	0.53	${0.94}_{-0.13}^{+0.41}$	${0.65}_{-0.13}^{+0.41}$	-1.2${}_{-0.1}^{+0.2}$	1.46 ± 0.08
A1664	0.128	113.8d	58.06	${100}_{-10}^{+50}$	${100}_{-10}^{+50}$	${1.00}_{-0.06}^{+0.16}$	2.21 ± 0.04
A1689	0.184	-0.28j	−0.28	<0.28	<0.28	<-1.50	2.31 ± 0.1
A1795	0.063	1.13b	10.63	${19}_{-3}^{+8}$	${19}_{-3}^{+8}$	${0.25}_{-0.06}^{+0.16}$	2.27 ± 0.02
A1991	0.059	4a	4.10	${7.3}_{-1.0}^{+3.2}$	${7}_{-1.0}^{+3.2}$	−0.17${}_{-0.06}^{+0.16}$	1.58 ± 0.04
A2052	0.035	1.8a	1.85	${3.3}_{-0.4}^{+1.4}$	${3.0}_{-0.4}^{+1.4}$	−0.55${}_{-0.07}^{+0.17}$	1.66 ± 0.02
A2142	0.090	0.02b	0.40	${0.72}_{-0.10}^{+0.31}$	${0.39}_{-0.10}^{+0.31}$	−1.4${}_{-0.1}^{+0.3}$	1.73 ± 0.07
A2199	0.030	2.7d	1.38	${2.5}_{-0.3}^{+1.1}$	${2.4}_{-0.3}^{+1.1}$	−0.65${}_{-0.06}^{+0.16}$	1.68 ± 0.05
A2204	0.152	159.4d	81.33	${150}_{-20}^{+60}$	${150}_{-20}^{+60}$	${1.1}_{-0.1}^{+0.2}$	2.7 ± 0.01
A2244	0.097	0.333i	0.33	${0.59}_{-0.08}^{+0.26}$	${0.45}_{-0.08}^{+0.26}$	-1.4${}_{-0.1}^{+0.2}$	1.47 ± 0.02
A2261	0.224	1.318i	1.32	${2.4}_{-0.3}^{+1.0}$	${2.0}_{-0.3}^{+1.0}$	-0.72${}_{-0.07}^{+0.18}$	2.1 ± 0.1
A2390	0.230	109a	111.00	${200}_{-30}^{+90}$	${200}_{-30}^{+90}$	${1.3}_{-0.1}^{+0.2}$	2.18 ± 0.06
A2597	0.085	29.7j	53.84	${96}_{-13}^{+42}$	${96}_{-13}^{+42}$	${0.96}_{-0.06}^{+0.16}$	2.49 ± 0.05
A2626	0.057	1.1d	0.56	${1.0}_{-0.1}^{+0.4}$	${0.90}_{-0.13}^{+0.44}$	−1.1${}_{-0.1}^{+0.2}$	1.21 ± 0.06
A3112	0.072	7.1a	7.28	${13}_{-2}^{+6}$	${12}_{-2}^{+6}$	${0.074}_{-0.065}^{+0.163}$	1.93 ± 0.05
A3581	0.022	2.4a	2.47	${4.4}_{-0.6}^{+1.9}$	${4.3}_{-0.6}^{+1.9}$	−0.39${}_{-0.06}^{+0.16}$	1.35 ± 0.22
A4059	0.048	4.1a	4.21	${7.5}_{-1.0}^{+3.3}$	${7.0}_{-1.0}^{+3.3}$	−0.18${}_{-0.07}^{+0.17}$	1.09 ± 0.06
Cygnus A	0.056	28.4j	21.32	${38}_{-5}^{+17}$	${38}_{-5}^{+17}$	${0.56}_{-0.06}^{+0.16}$	2.15 ± 0.01
Hercules A	0.154	29a	29.63	${53}_{-7}^{+23}$	${52}_{-7}^{+23}$	${0.70}_{-0.06}^{+0.16}$	1.75 ± 0.01
Hydra A	0.055	13a	13.34	${24}_{-3}^{+10}$	${24}_{-3}^{+10}$	${0.35}_{-0.06}^{+0.16}$	2.04 ± 0.02
MKW3S	0.045	−14.7f	0.95	${1.7}_{-0.2}^{+0.7}$	${1.6}_{-0.2}^{+0.7}$	−0.80${}_{-0.06}^{+0.16}$	1.36 ± 0.05
MS 0735.6+7421	0.216	9g	93.40	${170}_{-20}^{+70}$	${170}_{-20}^{+70}$	${1.2}_{-0.1}^{+0.2}$	2.42 ± 0.08
MS 1455.0+2232	0.259	29g	453.87	${810}_{-110}^{+350}$	${810}_{-110}^{+350}$	${1.9}_{-0.1}^{+0.2}$	2.78 ± 0.02
NGC 4782	0.013	1.02c	0.78	${1.4}_{-0.2}^{+0.6}$	${1.1}_{-0.2}^{+0.6}$	−0.99${}_{-0.08}^{+0.19}$	0.23 ± 0.03
NGC 5044	0.009	0.54a	0.56	${0.99}_{-0.13}^{+0.43}$	${0.85}_{-0.13}^{+0.43}$	−1.1${}_{-0.1}^{+0.2}$	1.94 ± 0.03
PKS 0745-191	0.103	63a	64.51	${120}_{-20}^{+50}$	${110}_{-20}^{+50}$	${1.0}_{-0.1}^{+0.2}$	2.89 ± 0.01
RXC J0352.9+1941	0.110	62a	63.47	${110}_{-20}^{+50}$	${110}_{-20}^{+50}$	${1.0}_{-0.1}^{+0.2}$	2.3 ± 0.03
RXC J1459.4-1811	0.230	240a	244.39	${440}_{-60}^{+190}$	${440}_{-60}^{+190}$	${1.6}_{-0.1}^{+0.2}$	2.48 ± 0.04
RXC J1524.2-3154	0.100	46a	47.11	${84}_{-11}^{+37}$	${84}_{-11}^{+37}$	${0.90}_{-0.06}^{+0.16}$	2.23 ± 0.01
RXC J1558.3-1410	0.100	22a	22.53	${40}_{-5}^{+17}$	${40}_{-5}^{+17}$	${0.58}_{-0.06}^{+0.16}$	2.1 ± 0.03
RXJ0439+0520	0.208	62a	63.20	${110}_{-10}^{+50}$	${110}_{-10}^{+50}$	${1.0}_{-0.1}^{+0.2}$	2.39 ± 0.23
RXJ1504.1-0248	0.215	475.595i	475.60	${850}_{-110}^{+370}$	${850}_{-110}^{+370}$	${1.9}_{-0.1}^{+0.2}$	3.29 ± 0.08
RXJ1539.5-8335	0.073	30a	30.76	${55}_{-7}^{+24}$	${55}_{-7}^{+24}$	${0.72}_{-0.06}^{+0.16}$	2.19 ± 0.05
RXJ1720.1+2638	0.164	12.7d	6.48	${12}_{-2}^{+5}$	${11}_{-2}^{+5}$	${0.038}_{-0.063}^{+0.159}$	2.63 ± 0.03
RXJ2129.6+0005	0.235	32a	32.58	${58}_{-8}^{+25}$	${57}_{-8}^{+25}$	${0.74}_{-0.06}^{+0.16}$	2.36 ± 0.33
Srsic 159-03	0.058	1.06b	8.53	${15}_{-2}^{+7}$	${15}_{-2}^{+7}$	${0.15}_{-0.06}^{+0.16}$	2.37 ± 0.02
Zw 2701	0.210	8.7d	4.44	${7.9}_{-1.1}^{+3.4}$	${7.8}_{-1.1}^{+3.4}$	−0.13${}_{-0.06}^{+0.16}$	1.81 ± 0.27
Zw 3146	0.290	704.7d	359.55	${640}_{-90}^{+280}$	${640}_{-90}^{+280}$	${1.8}_{-0.1}^{+0.2}$	2.87 ± 0.11
Perseus	0.018					1.85 ± 0.28k	2.67 ± 0.05

Notes. Column 1: system name. Column 2: redshift. Column 3: literature Hα measurements, with varying measurements quoted (e.g., Hα flux versus luminosity) and differing cosmologies (values of H₀, Ω_Λ, and Ω_m ) before homogenization to a Hα luminosity (i.e., Column 4): ^a Hamer et al. (2016), ^b McDonald et al. (2010), ^c Lakhchaura et al. (2018), ^d Crawford et al. (1999), ^e Gomes et al. (2016), ^f Buttiglione et al. (2009), ^g Donahue et al. (1992), ^h Heckman et al. (1989), ⁱ Wang et al. (2010), ^j Cavagnolo et al. (2009), and ^k Mittal et al. (2015). Column 5: Hα luminosity after intrinsic and Galactic extinction correction. Column 6: Hα luminosity after correcting for evolved stellar population contribution (planetary nebulae, supernova remnants, AGB and HB stars, etc.). Column 7: Hα SFR. Column 8: maximal (i.e., "classical") ICM cooling rate, from McDonald et al. (2018).

Download table as:  ASCIITypeset image

Both SFR and ${\dot{M}}_{\mathrm{cool}}$ for Figure 3 range over several orders of magnitude, and we see that our sample of starbursting BCGs lies at the extreme ends of both SFR and ${\dot{M}}_{\mathrm{cool}}$, with much higher cooling efficiencies (defined as the ratio ${\epsilon }_{\mathrm{cool}}\,=\,\mathrm{SFR}/{\dot{M}}_{\mathrm{cool}}$) than the reference systems. For the reference systems, we find a combined average cooling efficiency of 1.3 ± 0.1%, with a log-normal scatter of 0.65 dex, demonstrating the well-known 2 orders of magnitude suppression of the cooling flow problem. In contrast, our sample of extreme cooling clusters (N = 7) has an average cooling efficiency of 20% ± 13%. When binning the data points by cooling rate, one can readily see that systems in higher cooling rate bins have higher average cooling efficiencies, as shown by the colored histograms corresponding to the same colored data points in the scatterplot in Figure 3. Motivated by this trend, we fit the SFR versus ${\dot{M}}_{\mathrm{cool}}$ plot with a total least-squares regression using the python software package LinMix (Kelly 2007). LinMix ¹⁶ uses a hierarchical Bayesian approach to linear regression for data with both x- and y-errors, as well as robustly accounting for censored data (i.e., upper limits). We find a relation between cooling and star formation of $\mathrm{log}(\mathrm{SFR})\,=\,(1.66\,\pm \,0.17)\,\mathrm{log}({\dot{{\rm{M}}}}_{\mathrm{cool}})\,+(-3.22\pm 0.38)$ with an even lower intrinsic scatter of 0.39 ± 0.09 dex compared to the 0.65 dex from averaging over the entire sample and assuming a constant cooling efficiency. This relation tells us that at higher cooling rates, star formation proceeds with greater efficiency, consistent with McDonald et al. (2018) and Fogarty et al. (2015), though the latter found a steeper relation.

The steeper slope found by Fogarty et al. (2015) may be attributed to a few different factors. Figure 3 includes a large number of cool cores with ${\dot{M}}_{\mathrm{cool}}\lt 100$ M_⊙yr⁻¹ and a measured SFR, whereas Fogarty et al. (2015) does not. Additionally, the Fogarty et al. (2015) relation was based on a slightly different cooling rate definition than ours, where they measured ${\dot{M}}_{\mathrm{cool}}$ within a fixed 35 kpc aperture as well as one where t_cool/t_ff < 20. Furthermore, the CLASH sample considered in Fogarty et al. (2015) had a mean redshift of 〈z〉 = 0.392 and a complicated selection function, while the linear fit in Figure 3 was performed over many more systems, spanning from 0 ≲ z ≲ 0.5 (∼5 Gyr in lookback time), with an average redshift of 〈z〉 = 0.183. For comparison, we perform a BCES orthogonal fit to our data and find a best-fit relation of $\mathrm{log}(\mathrm{SFR})=(2.08\pm 0.15)\,\mathrm{log}({\dot{{\rm{M}}}}_{\mathrm{cool}})+(-4.12\pm 0.38)$, in closer agreement with the slope of Fogarty et al. (2015). We note that this steeper slope does not affect the interpretation of our results in the discussion below (Section 4). We see in the right panel of Figure 3 a slight redshift trend in the cooling efficiency histograms, where higher-redshift systems on average have higher _cool. While not highly significant, this redshift trend could likely be due to observational bias, where it is harder to find less-massive systems at higher redshifts. Alternatively, if real, this trend could be due to a transition between quasar-mode feedback generally observed more in higher-z systems, to radio-mode feedback that is often characteristic of more low-redshift systems on average (e.g., Churazov et al. 2005; Russell et al. 2013; Sadowski & Gaspari 2017). Our current sample is not suited for weighing in on this trend, but we will investigate in a future paper whether there is any redshift evolution of _cool in a larger, more complete, and unbiased Sunyaev–Zel'dovich (SZ)–selected sample (M. S. Calzadilla et al. 2022, in preparation).

The steeper-than-unity relation between SFR and ICM cooling rate presented in Section 3.2 implies that multiphase condensation gradually becomes more prevalent as the maximal cooling rate increases. In this section, we attempt to explain only why the conversion from hot (10⁷ K) to warm (10⁴ K) phases becomes more efficient with cooling rate, and not whether the conversion from warm gas to stars in these systems is more efficient. In cool-core clusters, the ICM density peaks sharply at the cluster center, where the condensing material accumulates within the BCG. This condensing material should eventually form cold molecular gas reservoirs that fuel star formation and accrete onto the central SMBHs. A large amount of the cooling in the most extreme cooling systems may be nonradiative (e.g., mixing/dust cooling). Regardless, despite clear evidence of feedback from the AGN (e.g., X-ray cavities), these outbursts are unable to suppress cooling to the same degree as in systems with lower cooling rates. One way to understand why this is the case is to contrast the growth rates of the SMBHs versus that of the cluster halos they reside in. Using various empirical scaling relations, McDonald et al. (2018) demonstrated that the accretion rate of SMBHs (${\dot{M}}_{\bullet }$) scaled by the Eddington rate is related to the ICM cooling rate as ${\dot{M}}_{\bullet }/{\dot{M}}_{\mathrm{Edd}}\propto {\dot{M}}_{\mathrm{cool}}^{1.87}$, which is consistent with accretion rate data from Russell et al. (2013), and is supported by the tight positive scalings between the SMBH mass and hot halo properties (Gaspari et al. 2019). This correlation implies that more massive clusters, with very high cooling rates (${\dot{M}}_{\mathrm{cool}}\gtrsim {10}^{3}$ M_⊙yr⁻¹), should have a central AGN accreting closer to the Eddington rate than in low-mass systems. However, while the left-hand side of the above relationship can "saturate as the SMBH growth rate is capped by the Eddington limit, the right-hand side has no analogous constraint. Galaxy cluster halos grow via mergers and accretion of smaller halos, which can relatively quickly increase the available amount of cooling material. In the most massive galaxy clusters, where AGN accretion approaches the Eddington rate, we might then expect disproportionately undermassive SMBHs and for ICM cooling to increasingly outpace the feedback at higher cooling rates, thus steepening the SFR–${\dot{M}}_{\mathrm{cool}}$ relation. This is a testable hypothesis, though the direct measurement of SMBH masses and accretion rates in BCGs, especially those with quasar-hosting systems, is difficult (e.g., McConnell & Ma 2013).

Another consequence of the higher radiative efficiency of an AGN in more strongly cooling halos is the resulting dominant mode of AGN feedback, which has an impact on how well the AGN energy can couple to the cooling ICM. As radiative efficiency ${\dot{M}}_{\bullet }/{\dot{M}}_{\mathrm{Edd}}\to 1$, a higher fraction of the AGN's accretion energy gets released in the form of radiation rather than in mechanical outflows via jets (Churazov et al. 2005; Russell et al. 2013; Sadowski & Gaspari 2017). This transition from mechanical to radiative feedback is gradual and incomplete, meaning that it is not a sudden switch where radio jets turn off. Phoenix, H1821, and IRAS09104 are excellent examples of quasar-hosting systems that also have jet-inflated bubbles. Observationally, this radiative energy output begins to dominate gradually as the black hole's accretion rate approaches and exceeds a few percent of the Eddington rate, i.e., ${\dot{M}}_{\bullet }\gtrsim 0.1{\dot{M}}_{\mathrm{Edd}}$, corresponding to a cooling rate of ∼1000 M_⊙yr⁻¹ (see Figure 8 in McDonald et al. 2018). Above this accretion rate is where the quasar-hosting systems in our extreme cooling sample reside (see Figure 12 in Russell et al. 2013). We argue that it is no coincidence that the systems in our extreme cooling sample, particularly the quasar-hosting clusters, are also among the most highly efficient at converting hot gas into stars. An AGN in the radiative mode may allow more cooling to occur since the hot atmosphere it is embedded in is optically thin to radiation, making it less capable of coupling large amounts of its accretion energy to its surroundings and quenching cooling compared to a mechanical mode AGN. By contrast, radiatively inefficient, mechanical mode AGN can heat their surroundings via a number of simultaneous channels since their far-reaching jets can inflate bubbles, which do work when expanding against their surroundings, as well as create shocks and sound waves, release cosmic rays, and mechanically increase the turbulence in the ICM (e.g., Soker et al. 2001; Churazov et al. 2001; Reynolds et al. 2002; Zhuravleva et al. 2014; Gaspari 2015; Yang & Reynolds 2016; Hillel & Soker 2017; Li et al. 2017; Yang et al. 2019). The fact that jets are "on" for a large fraction of the time (∼60%–70%; e.g., Dunn & Fabian 2006; Bîrzan et al. 2012) is further evidence that specifically mechanical feedback is needed to regulate star formation and prevent runaway cooling. It would additionally suggest that radiatively efficient feedback is at least predominantly operating in, if not responsible for, those systems where cooling is overpowering heating.

The Eddington limit might play an even more critical role in systems with high-mass cores considering variations in the accretion flow onto the SMBH. If the average accretion rate $\langle {\dot{M}}_{\bullet }\rangle \sim 0.1{\dot{M}}_{\mathrm{Edd}}$, and accretion is clumpy and chaotic rather than steady, then any small clump of material that gets accreted will abruptly spike the accretion rate to the Eddington limit and suppress the most energetic outbursts via radiation pressure. In other words, at high average accretion rates, scatter can no longer be log normal because one side is truncated by the Eddington limit, while the other side is not. Thus, the average effect of feedback may be reduced. Moreover, the role of additional Compton cooling induced by quasars in the central few kiloparsecs is to suppress feedback effectiveness further, but it should only be an appreciable effect for an unobscured quasar like H1821 (Russell et al. 2010; McDonald et al. 2019). Given this potential time dependence, it is natural to ask whether the extreme cooling we see in our sample is a short-lived phenomenon that occurs in most cool-core clusters or if it is characteristic of only a small subset of cool-core clusters. Somboonpanyakul et al. (2021) searched for Phoenix-like systems misidentified in the ROSAT survey (Voges et al. 1999) as X-ray bright point sources, and concluded that similar starbursting BCGs have an occurrence rate of ≲1%. Prior to that, McDonald et al. (2019) used deep Chandra data to show that the Phoenix cluster is the strongest example of a potential runaway cooling flow out of ∼180 cool-core clusters ranging over 9 Gyr in time. Such an occurrence rate implies that if extreme cooling as seen in Phoenix is a common phenomenon, then it must only have a duration of roughly ∼50 Myr. This is consistent with simulations (e.g., Prasad et al. 2015, 2020) showing that episodes of high cooling rates and SFRs should last for ≲100 Myr in any cool-core cluster. However, the idea of cool-core cycles may be inconsistent with our findings of a slope steeper than unity in the SFR–${\dot{M}}_{\mathrm{cool}}$ plot as well as the decreasing scatter with higher cooling rates shown in Figure 3. For instance, the lack of clusters with both ${\dot{M}}_{\mathrm{cool}}$ ≳ 1000 M_⊙yr⁻¹ and SFR ≲30 M_⊙yr⁻¹, which should be relatively easy to detect and are missing even in the more complete sample compiled by McDonald et al. (2018), tells us that our extreme cooling sample is probably not a sample of clusters caught at an opportune time (i.e., at a cooling extremum in the cool-core cycles described in Prasad et al. 2015, 2020). Still, the issue of how these extreme systems then stop cooling and quench star formation requires further investigation. It may be that with the more efficient conversion from cooling ICM to star formation that happens at high cooling rates (${\dot{M}}_{\mathrm{cool}}\,\gtrsim \,1000$ M_⊙yr⁻¹), and thus closer to Eddington accretion rates (${\dot{M}}_{\bullet }\gtrsim 0.1{\dot{M}}_{\mathrm{Edd}}$), that the initially undersized SMBH eventually grows enough to increase the Eddington accretion rate itself, making the accretion fall into the sub-Eddington, mechanical feedback-dominated regime again. Indeed we do see ultramassive black holes up to several 10¹⁰ M_⊙, which may push down the Eddington rates (see Gaspari et al. 2019). If we assume a CCA evolution, we often see spikes in accretion up to the Eddington regime lasting a few megayears, but given the chaotic nature of the feedback, the grown SMBH quickly self-regulates down to lower rates with a flickering time power spectrum described as in Gaspari et al. (2017). More observational studies or simulations testing this hypothesis could be a promising path forward.

The onset of thermal instabilities in the hot ICM is currently a contentious issue. Heating via mechanical feedback largely suppresses cooling in cool-core clusters, but some cooling still happens, and moreover it is aided by this same self-sustaining feedback loop. Using the maps in Figure 2, we have seen where this cooling happens. Now, we can use these maps to compare the extent of the [OII]-emitting gas to recently established metrics that attempt to explain the details of how thermal instabilities develop in the presence of AGN feedback. In this section, we will explore whether the measured extent of cool nebular gas in these maps correlates with X-ray derived radii that indicate cooling instability. Extent measurements like these can be affected by projection along the line of sight, as well as observing depth, both of which lead to these measurements being lower limits of a true multiphase cooling threshold radius. To counteract these limitations, we add to the seven systems in our extreme cooling sample those of McDonald et al. (2010) for which we obtained Hα maps. We also utilize the Hα map of Perseus (NGC1275; Conselice et al. 2001), one of the closest, brightest, and most well-studied clusters. The network of Hα filaments in Perseus may be representative of the range of extents and angles to our line of sight that we could expect in other more distant systems. One caveat in Perseus (and possibly others) is that not all of the Hα emission is associated with star formation (e.g., Canning et al. 2010, 2014), with collisional heating being a potential ionization source instead (Ferland et al. 2009). Even so, these Hα and [OII] maps of Perseus and others still trace where we see warm ∼10⁴ K multiphase gas that has cooled out of an unstable hot atmosphere.

Rather than measure a single maximum extent of nebular gas in all of these systems, we measure the maximum extent in 10 sectors, evenly spaced in azimuth and centered on the BCG position to find an average extent, 〈R〉_Hα,[OII]. This more robust measurement further reduces the sensitivity that a single maximum extent has to small noise blobs at large radii. Additionally, it is more fair to compare an azimuthally averaged extent to the azimuthally averaged t_cool and t_cool/t_ff profiles we will analyze in Section 4.2.2. These extent measurements are all listed in Table 5, and plotted against different X-ray diagnostics of thermal instability in Figure 4. In each of these panels, the shared y-axis measurements of 〈R〉_Hα,[OII] have values and uncertainties (in solid colors) determined from the median and interquartile range (i.e., 25th and 75th) of the extent distributions, respectively. The larger, dashed, gray y–error bars on each of these data points depict the minimum and maximum extents across all azimuthal sectors for each cluster, which serves to encode the asymmetry of certain systems like A1795, for instance, which are extended along mostly one direction or axis. In the following, we examine the relationship between these average extents and various indicators of ICM instability.

Zoom In Zoom Out Reset image size

Table 5. Various Extent Data Used in Figure 4

Name	rbubble		r @ (tcool/tff)inflection	r @ tcool = 1 Gyr	RHα,[OII] (kpc)
	(center) (kpc)	(center+radius) (kpc)	(kpc)	(kpc)	(median)	(min., max.)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
Phoenix	36.3	50.0	<3.3 a	${59}_{-1}^{+2}$	${27}_{-7}^{+17}$	(15, 63)
H1821	24.2	33.3	<30 a	${38}_{-3}^{+2}$	${31}_{-3}^{+5}$	(24, 46)
IRAS09104	41.7	63.3	${27}_{-8}^{+5}$	${45}_{-2}^{+7}$	${23}_{-3}^{+4}$	(14, 31)
A1835	18.8	25.0	${50}_{-6}^{+10}$	${37}_{-2}^{+1}$	${19}_{-4}^{+2}$	(10, 24)
MACS1931	25.3	35.4	${55}_{-7}^{+7}$	46 ± 1	29 ± 3	(23, 36)
RXJ1532	33.8	50.2	${60}_{-17}^{+16}$	${49}_{-3}^{+7}$	${25}_{-5}^{+4}$	(15, 35)
RBS797	32.3	49.2	${54}_{-3}^{+1}$	52 ± 1	24 ± 3	(18, 33)
A85	21.3	28.8	21 ± 6	17 ± 3	${1.7}_{-0.4}^{+0.6}$	(1.2, 3.2)
A133	32.4	61.5	26 ± 3	18 ± 1	${2.4}_{-0.7}^{+0.2}$	(1.2, 4.7)
A478	9.0	13.4	30 ± 6	27 ± 4	${7.7}_{-1.2}^{+0.7}$	(5.3, 11)
A496			21 ± 4	17 ± 2	${4.0}_{-0.6}^{+1.3}$	(3.0, 7.4)
A1644			23 ± 3	7 ± 1	${6.2}_{-4.1}^{+5.1}$	(0.9, 17)
A1795	18.5	30.0	64 ± 13	20 ± 5	${8.8}_{-2.2}^{+2.4}$	(4.9, 55)
A2052	11.2	20.4	15 ± 3	21 ± 1	${6.3}_{-1.4}^{+2.8}$	(3, 16)
A2597	22.6	31.1	36 ± 6	28 ± 3	${14}_{-3}^{+7}$	(8, 27)
A4059	22.7	37.1	26 ± 8	13 ± 9	${3.2}_{-0.5}^{+0.7}$	(2.1, 5.1)
Srsic 159-03	26.4	45.3			${9.0}_{-5.5}^{+5.8}$	(1.5, 32)
Perseus b				37 ± 1	${27}_{-5}^{+4}$	(13, 63)

Notes. Column 1: system name. Columns 2 and 3: X-ray cavity/bubble data measured from our archival Chandra X-ray images (systems above horizontal rule), or taken from Diehl et al. (2008). Distance is measured to the bubble center (in Column 2) as well as to the leading edge (center+radius; Column 3) of the bubble. Column 4: inflection point in the t_cool/t_ff profile, modeled as a double power law with a floor (i.e., inner power-law slope of zero). Column 5: radius where cooling time falls below 1 Gyr. Column 6: median and interquartile range (25th and 75th percentiles) of the distribution of optical filament extents measured from our [OII] maps (i.e., systems above horizontal rule; see Figure 2) or the Hα measurements of McDonald et al. (2010). Column 7: minimum and maximum optical filament extents.

^a The innermost radial bin from the X-ray data is quoted as an upper limit since the t_cool/t_ff profile is consistent with a single power law, with no resolved inflection point in the profile. ^b Perseus Hα data quoted from Conselice et al. (2001) and Fabian et al. (2003), and t_cool data from Dunn & Fabian (2006). No bubble data are quoted for this system due to the multiple X-ray cavity pairs seen in Perseus, making the association between a particular cavity and the filaments difficult.

Download table as:  ASCIITypeset image

4.2.1. 〈R〉_Hα,[OII] versus Bubbles

4.2.2. 〈R〉_Hα,[OII] versus Cooling and Freefall Time Profiles

In addition to the "stimulated feedback" model described above, other models predict that thermal instabilities can be produced by the cooling of the turbulent ICM when and where it becomes multiphase, i.e., where the specific entropy of the gas falls below some threshold value (e.g., Cavagnolo et al. 2008), resulting in "precipitation" of low angular momentum gas clouds onto the central SMBH. Precipitation models (e.g., Gaspari et al. 2012; Voit et al. 2015) predict that these thermal instabilities develop within the transition radius between an outer baseline ICM entropy profile due to cosmological structure formation and an inner profile induced by CCA and feedback.

While models predict that thermally unstable cooling happens where t_cool/t_ff < 10, most t_cool/t_ff profiles of cool-core clusters, with the exception of Phoenix, do not fall significantly below this threshold, as seen in Hogan et al. (2017). Instead, Hogan et al. (2017) argued that because mass profiles within the typical radii where t_cool/t_ff approaches a minimum can be approximated by an isothermal sphere, and that entropy profiles within these typical radii are consistent with a single power law, then t_cool/t_ff profiles should flatten out toward smaller radii. Furthermore, in most of the 33 Hα emitting clusters studied by Hogan et al. (2017), the min(t_cool/t_ff) values were measured from an annulus just outside of a single inner core annulus with a noisy t_cool/t_ff measurement, showing that these measurements are typically not well resolved. For these reasons, we modeled each of the t_cool/t_ff profiles for our extreme cooling sample (see Section 2.3.1) as well as those from the Hogan et al. (2017) clusters that have corresponding Hα extent data from McDonald et al. (2010) with the assumption of a floor value rather than a minimum. To do so, we fit two power laws to the t_cool/t_ff profiles, fixing the slope of just the innermost power law to zero and allowing both normalizations to vary. We performed each of these fits over 1000 bootstrap iterations to obtain a distribution of measurements of where the two power laws cross, and plot these inflection radii ($R{({t}_{{\rm{c}}}/{t}_{\mathrm{ff}})}_{\mathrm{inflection}}$) against Hα and [OII] extents in Figure 4 (middle panel). We again find a weak correlation (ρ = 0.45, p_ρ = 0.07 ; r = 0.31, p_r = 0.22), but note that the correlation strength is driven down largely by two upper limits on the inflection radius. These two upper limits include Phoenix, whose t_cool/t_ff profile has no discernible floor or minimum in the Chandra data, and H1821, whose bright point source in the X-ray observations prevented a resolved measurement of $R{({t}_{{\rm{c}}}/{t}_{\mathrm{ff}})}_{\mathrm{inflection}}$. The fact that ρ > r could be an indication of a slightly nonlinear relationship between ($R{({t}_{{\rm{c}}}/{t}_{\mathrm{ff}})}_{\mathrm{inflection}}$) and extent of multiphase gas. Just as in the bubble correlation plot discussed above, we see that most systems have an inflection point in their t_cool/t_ff profile that surrounds the average volume within which we see cooling filaments.

Some studies have found that the scatter of t_ff values is significantly smaller than the range of t_cool values at either a fixed radius of 10 kpc (Hogan et al. 2017) or at altitudes where t_cool/t_ff is at a minimum (McNamara et al. 2016). These findings suggest that the local gravitational effects encoded in t_ff do not add any predictive power for the onset of thermal instabilities over t_cool alone (Hogan et al. 2017). Some difficulties also arise from calculating t_ff profiles. In light of this, we also compare in Figure 4 (right panel) where the individual modeled t_cool profiles of Hogan et al. (2017) first decrease past 1 Gyr versus the average multiphase cooling radius. Hogan et al. (2017) showed that the deprojected t_cool profiles of clusters with no nebular emission do not go below this 1 Gyr threshold at a radius of 10 kpc, with the exception of A2029. We show that there is a very strong correlation between t_cool and the average extent of cooling (ρ = 0.88, p_ρ = 1.9 ×10⁻⁶; r = 0.89, p_r = 8.1 × 10⁻⁷). The addition of our seven new extent measurements (plotted in blue) is especially helpful in establishing this correlation. Once more, we see that all of the average extents in the right panel of Figure 4 lie below the one-to-one line, indicating that we observe multiphase gas out to radii where the cooling time is ≲1 Gyr. It could be argued that this cooling time threshold is somewhat arbitrary, as a much shorter threshold (e.g., 0.1 Gyr) would move all of the data points to the left, possibly above the one-to-one line in Figure 4. This threshold should vary over the mass ranges of rich clusters down to poor groups, where central cooling times can be up to 10× shorter, highlighting the importance of some mixing timescale for normalization (e.g., t_ff or t_eddy). Despite this, it is worth noting that the strong correlation between 〈R〉_Hα,[OII] and R(t_cool = 1 Gyr) should persist and always arise in any scenario in which multiphase gas is supplied locally by the cooling of hot ambient gas, which nicely supports the notion that the presence of multiphase gas is linked to cooling of ambient gas, regardless of the mechanism.

The strong correlation between 〈R〉_Hα,[OII] and R(t_cool = 1 Gyr) in the right panel of Figure 4 is particularly interesting as it hints at the best estimate yet of where thermally unstable cooling ensues. By trying a range of cooling time thresholds below 1 Gyr, we find that measuring the radius at which t_cool = 0.5 Gyr brings the data points closest to a one-to-one correspondence with the extent of optical filaments, as shown in Figure 5. This tight relation suggests that, on average, when the ICM reaches a cooling time of 500 Myr or shorter, thermal instabilities will develop, leading to extended filaments of multiphase gas.

Zoom In Zoom Out Reset image size

More broadly, Figure 4 shows for the first time that multiphase gas on average resides beneath the altitudes where we see X-ray cavities, the precipitation limit marked by a change in slope in t_cool/t_ff profiles, and where the average cooling time is shorter than 1 Gyr. No matter which indicator one chooses, they all circumscribe the average volume within which the ICM is thermally unstable to multiphase cooling. The fact that the average extents of filaments do not lie on the one-to-one line with the altitudes defined by the various X-ray measurements is unsurprising as the snapshot in time at which we observe filaments at optical wavelengths could be drastically different from the onset of when the local cooling cascade in the ICM begins. In hydrodynamic simulations, cold gas is seen to quickly fall out of pressure equilibrium with hot gas (e.g., Qiu et al. 2021). Turbulence could also play a large role not only in the altitude to which we observe multiphase gas, but also in the diverse morphology of [OII] nebulae that we see in Figure 2 for instance. It could be that measuring the eddy timescale of the turbulent ICM (e.g., Gaspari et al. 2018) is more strongly correlated with the maximum extent of multiphase gas, but these measurements are difficult to make in practice, requiring filament kinematics. We may be able to more closely connect the hot (> 10⁷ K) gas to the intermediate warm (10^5.5 K) gas, and better assess how ICM cooling flows fuel star formation, using future observations of coronal emission lines with instruments like MIRI aboard the James Webb Space Telescope.