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

We present new, deep, narrow- and broad-band Hubble Space Telescope observations of seven of the most star-forming brightest cluster galaxies (BCGs). Continuum-subtracted [O II] maps reveal the detailed, complex structure of warm ($T \sim 10^4$ K) ionized gas filaments in these BCGs, allowing us to measure spatially-resolved star formation rates (SFRs) of ~60-600 Msun/yr. We compare the SFRs in these systems and others from the literature to their intracluster medium (ICM) cooling rates (dM/dt), measured from archival Chandra X-ray data, finding a best-fit relation of log(SFR) = (1.67+/-0.17) log(dM/dt) + (-3.25+/-0.38) with an intrinsic scatter of 0.39+/-0.09 dex. This steeper-than-unity slope implies an increasingly efficient conversion of hot ($T \sim 10^7$ K) gas into young stars with increasing dM/dt, or conversely a gradual decrease in the effectiveness of AGN feedback in the strongest cool cores. We also seek to understand the physical extent of these multiphase filaments that we observe in cluster cores. We show, for the first time, that the average extent of the multiphase gas is always smaller than the radii at which the cooling time reaches 1 Gyr, the tcool/tff profile flattens, and that X-ray cavities are observed. This implies a close connection between the multiphase filaments, the thermodynamics of the cooling core, and the dynamics of X-ray bubbles. Interestingly, we find a one-to-one correlation between the average extent of cool multiphase filaments and the radius at which the cooling time reaches 0.5 Gyr, which may be indicative of a universal condensation timescale in cluster cores.

Most of the baryonic matter in galaxy clusters is in the form of a hot (∼10 7 K), X-ray emitting plasma that permeates the space between the galaxies, known as the intracluster medium (ICM), leaving only a few percent of the mass budget to be found in stars. In so-called "cool core" clusters, as particles in the ICM interact and radiate away energy, they plunge deeper into the dark matter potential well of the cluster, eventually cooling out of the hot plasma phase until becoming cold enough to form stars. However, early studies of these multiphase cooling flows (e.g. Fabian 1994) found that only 1 − 10% of this cooling gas is actually observed to form stars. Apparently, most of the hot plasma in clusters is being kept hot by some heating source. This dilemma is referred to as the "cooling flow problem." In the past decade, feedback from active galactic nuclei (AGN) has emerged as the most likely heating source energetically capable of preventing the runaway cooling of the ICM (see reviews by McNamara & Nulsen 2012;Gaspari et al. 2020;Donahue & Voit 2022). In this scenario, the activity of an AGN is driven by the accretion of infalling material onto a supermassive black hole (SMBH), which then self-regulates its fuel supply via either radiation pressure (i.e. "quasar mode" feedback), or mechanical energy from an outburst that launches relativistic jets of plasma on tens of kpc scales ("radio mode" feedback). The tight correlation between the cooling luminosity of the ICM in relaxed clusters and the radio power (e.g. Pasini et al. 2021) as well as the outburst power -as measured by the work done by bubbles inflated by these relativistic jets as they expand against the ICM -establishes the SMBH as a thermostat that adds more heat to its enviroment as the central cooling time decreases, maintaining a gentle equilibrium (e.g. Hlavacek-Larrondo et al. 2015).
A wealth of observational evidence now backs up this AGN feedback model (see review by Fabian 2012). However, the discovery of the Phoenix cluster in 2012 has since prompted closer investigations into the limits of AGN feedback (McDonald et al. 2012b. This galaxy cluster hosts the most strongly starbursting brightest cluster galaxy (BCG) and the only known possible runaway cooling flow observed in any cool core galaxy cluster. In the most massive clusters, with classically-inferred maximal ICM cooling rates (Ṁ cool = M gas (r < r cool )/t cool ) on the order of >1000 M yr −1 , BCGs seem to be forming stars at rates of >10% of the cooling rate Molendi et al. 2016;Fogarty et al. 2015). In other words, these rare, rapidly-cooling systems host AGN that appear to be incapable of the roughly two orders of magnitude suppression of cooling that we find in most other systems, perhaps signalling a gradually progressing saturation of AGN feedback . Of particular interest in these extreme cooling systems and even in more modestly star-forming clusters is the role of environment. A central ICM entropy threshold of K 0 < 30 keV cm 2 , or where the cooling time is roughly t cool < 1 Gyr, has been shown to be a remarkably sensitive indicator of whether accretion onto the central SMBH fuels AGN activity or star formation "downstream" of the cooling flow (e.g. Nulsen 1986;Pizzolato & Soker 2005;Cavagnolo et al. 2008;Rafferty et al. 2008;Main et al. 2017;Hogan et al. 2017;Pulido et al. 2018). Now, in addition to the standard picture of AGN feedback where heating outbursts suppress the runaway cooling of the ICM, these same outbursts may enhance cooling via bubbles and turbulence, and allow the feedback loop to sustain itself.
Recent studies into how the ICM becomes thermally unstable to cooling and thus fuels the formation of stars and giant molecular clouds (e.g. Pulido et al. 2018) have focused on the importance of local gravitational acceleration. In particular, the ICM was originally thought to become unstable when the ratio of the cooling time t cool to the free-fall time t ff falls roughly below unity under the assumption of a static plane-parallel atmosphere (Cowie et al. 1980;Nulsen 1986). A more realistic three-dimensional atmosphere with turbulence should have a t cool /t ff threshold of about 10 (e.g. McCourt et al. 2012;Sharma et al. 2012;Gaspari et al. 2012). When and where this threshold is reached in the cluster atmosphere, that parcel of gas cannot persist without precipitating in a "rain" of cold clouds. This rain fuels feedback, which then heats up and decreases the density of the cluster core, thus increasing the cooling time (t cool ∝ n −1 e k B T 0.5 ) and t cool /t ff ratio. Eventually, this halts the precipitation which turns off the feedback and allows the atmosphere to eventually cool again (e.g. Li et al. 2015;. This selfregulating behavior has been reproduced successfully in various simulations, producing cluster atmospheres with 10 < min(t cool /t ff ) < 30 (e.g. Sharma et al. 2012;Gaspari et al. 2012Gaspari et al. , 2013Li & Bryan 2014;Li et al. 2015;Meece et al. 2017). Observations of cool core clusters with nebular emission, molecular gas, etc. have similarly been found to have inner t cool /t ff ratios approaching 10, with no examples significantly below this value (with the exception of the Phoenix cluster), suggesting that this value is not a threshold but rather a floor (e.g. Hogan et al. 2017;Pulido et al. 2018).
A closely related framework for instability in the ICM is that of "chaotic cold-accretion" (CCA), where merger or AGN-driven turbulence induces enhanced cooling (Gaspari et al. 2013Prasad et al. 2017;Wittor & Gaspari 2020). This turbulence seeds a population of cool clouds with a broad angular mo-mentum distribution. Clouds comprising the low-end tail of this distribution fall inwards towards the center of the cluster and fuel enhanced accretion onto the SMBH, relative to the approximately Bondi accretion found in more homogeneous media (Pizzolato & Soker 2005;Gaspari et al. 2013;Tabor & Binney 1993;Binney & Tabor 1995;Prasad et al. 2015). In this CCA model, the condition for thermally unstable cooling to occur is that the cooling time is small compared to the local dynamical or turbulent "eddy" time (t cool /t eddy 1; e.g. Gaspari et al. 2018), rather than an order of magnitude longer than the local freefall time (i.e. t cool /t ff 10). In the "stimulated feedback" model (e.g. Mc-Namara et al. 2016), thermally unstable cooling happens in-situ when warm (∼1 keV) X-ray gas is uplifted in the wake of AGN-inflated radio bubbles as they rise buoyantly out of the central cluster potential, which increases the infall time and promotes the condensation of cold gas (see also Pope et al. 2010). Evidence for this phenomenon has been provided in a number of systems where large reservoirs (∼10 10 M ) of cold (10 − 100 K) molecular gas, observed with the Atacama Large Millimeter Array (ALMA), is projected behind or draped around the location of X-ray cavities as seen by Chandra (e.g. A1835: McNamara et al. 2014, Phoenix: Russell et al. 2017a, MACS 1931.8-2634: Ciocan et al. 2021, A1664: Calzadilla et al. 2019, PKS 0745-191: Russell et al. 2016, A2597: Tremblay et al. 2018, A1795: Russell et al. 2017b, 2A 0335+096: Vantyghem et al. 2016). On the other hand, several systems show cold gas in hot halos, even without a direct correlation with bubbles (e.g. Temi et al. 2018;Olivares et al. 2019;Rose et al. 2019;Maccagni et al. 2021;North et al. 2021;McKinley et al. 2021).
To distinguish between the various models that explain the onset of nonlinear thermal instabilities in the ICM, one can look at where those thermal instabilities develop. For instance, stimulated feedback theory predicts that one ought to see condensation only up to a maximum altitude where radio bubbles could uplift lower entropy and lower altitude material, or compress higher-altitude material at the bubble's leading edge. Meanwhile, precipitation theory predicts that this condensation of cold gas should happen where t cool /t ff 10, or equivalently where the ICM entropy profile changes its slope due to heating from the AGN. One of the most successful ways to identify where multiphase cooling occurs has been via narrow-band surveys for Hα emission in low-redshift (z 0.1) clusters (e.g. Heckman et al. 1989;Crawford et al. 1999;Hatch et al. 2007;McDonald et al. 2010;Hamer et al. 2016). As we have discussed that the ICM becomes multiphase when the central cool-ing time falls below ∼1 Gyr (e.g. Cavagnolo et al. 2008), and that the cooling time of the ICM where it is coincident with the presence of Hα emission is ∼4× shorter than the surrounding ICM (e.g. McDonald et al. 2010), we can use the morphology of extended Hα nebulae to map out where the ICM is undergoing multiphase condensation.
While Hα emission traces warm (∼10 4 K) gas that has already cooled but not yet formed stars or accreted onto the central SMBH, few observatories are capable of probing the redshifted Hα wavelengths of more distant clusters (z 0.3), especially with the required angular resolution to study these thin ( 1 kpc) filamentary nebulae in detail. To address this gap in the literature, we present in this study ∼40 orbits of new Hubble Space Telescope (HST) data on seven of the most extremely star-forming ( 100 M yr −1 ) BCGs known (see e.g. Mc-Donald et al. 2018). Given the relatively high average redshift of these systems ( z avg = 0.38), we obtained deep, high-angular resolution maps of [O ii] λλ3726,3729 emission rather than Hα, providing a cleaner and higher signal-to-noise view of the thermally unstable regions of these extreme BCGs. Both Hα and [O ii] have been used extensively as star formation indicators, and in this study we will assume that both probe regions of active star formation. Our new [O ii] maps will allow us to correlate line emitting gas to the cooling ICM, star forming clumps, and radio jets/bubbles. In §2 we describe our HST observations, observing strategy, and data reduction in detail, as well as the archival Chandra data that we use to investigate any correspondences between the warm (∼10 4 K) ionized [O ii] gas and the spatial and spectral properties of the hot (∼10 7 K) ICM. The resulting maps and star formation rates (SFRs) associated with these new observations are presented in §3, with their implications for the limits of AGN feedback and the onset of thermal instabilities discussed in § 4. Finally, we conclude with the take-away points in § 5. Throughout this paper, we assume a flat ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7. All measurement errors are 1σ unless noted otherwise.

Sample Selection
Our "extreme cooling sample" originally consisted of six of the most star-forming (SFR 100 M yr −1 ) and strongly-cooling (Ṁ cool 1000 M yr −1 ) 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  (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 (Ṁ cool ∼1000 M yr −1 ), 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 §3.2.

Optical/HST
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) which were previously published , as well as 2 orbits each from the CLASH survey (Postman et al. 2012) for MACS1931 (GO12456) and RXJ1532 (GO12454). Our observing strategy was designed to tune the ramp filters to be centered on the [O ii] λλ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 ([O ii] , [O iii] λλ4959,5007, and Hα) to better model and subtract the underlying continuum from the narrowband observations (described in more detail below in §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 [O ii] emission.
For each system, we used STScI's DrizzlePac software package 2 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 sub-pixel 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 3 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.

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 T 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' (Y) and 'old' (O) template fluxes: This equation can be inverted to solve for the vector , 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 ): Finally, we can subtract this predicted narrowband flux from the observed narrowband flux to get a continuumsubtracted [O ii] -only flux. This SED fitting and continuum subtraction is done pixel-by-pixel to get [O ii] -only maps for each of the clusters in our sample. The procedure is illustrated in Figure 1.

X-ray/Chandra
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 [O ii] maps. A summary of the observations used for this analysis is presented in Table 2. All Chandra observations were  , 16135, 16545, 19581, 548 19582, 19583, 20630, 20631, 20634, 20635, 20636, 20797 H1821+643 9398, 9845, 9846 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 blanksky 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.

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 onsource spectrum was fit using a combined spectral model PHABS*(APEC + BREMSS) + APEC, where the first APEC component is to model thermal emission from the opti-

IRAS 09104+4109
10'' = 57 kpc Figure 1. Illustration of the pixel-by-pixel spectral energy distribution (SED) fitting procedure we use to produce emissionline only maps (see Figure 2). For each object, we use narrowband observations tuned to the redshifted wavelength of [O ii] λλ3726,3729, bracketed by images of broadband filters chosen to be free from strong emission lines usually associated with star formation or AGN emission (e.g. Hα and [Oiii]). The flux in a single pixel from the bracketing broadband images is then used to construct an SED and is fit by a linear combination of a 10 Myr starburst and 6 Gyr passive elliptical galaxy template, allowing their normalizations to vary (see top right and bottom right panels, highlighting two distinct pixel fits). The combined, predicted template flux at the wavelength of the narrowband filter is then subtracted from the narrowband image at that pixel, resulting in a continuum-subtracted image after the procedure is run for every common pixel in the aligned images.
cally 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): 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 0 , 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 or-der 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) ≡ n e n p dV = η × 4π × 10 14 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): which is a modified double-beta model with a cusp, rather than a flat core (defined by n 0 , r c , α, β), a steeper outer profile slope (defined by r s , γ, and ), and a cool core component (defined by n 02 , r c2 , and β 2 ). 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 B T n −2/3 e ), cooling time (t cool = 3 2 (ne+np)kBT nenpΛ(kBT,Z) ), and freefall time (t ff = 2r/g). We used the cooling function Λ(k B T, Z) from Sutherland & Dopita (1993), as parameterized by Tozzi & Norman (2001) (see also Parrish et al. 2009;Guo & Oh 2008). 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 (NFW) profile (Navarro et al. 1997) 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.

[O ii] Maps
Our new, reduced broad-and narrow-band HST observations can be found in Figure 2. The continuumsubtracted [O ii] maps in the middle panels reveal with remarkable detail the intricate morphology of the warm (T ∼10 4 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 starforming material at optical wavelengths in each of these strongly cooling clusters to date (with the exception of  Fabian et al. 2003), with R Hα ≈ 40 kpc.
The morphologies of the [O ii] maps appear to be extending either in various directions (as is the case in Phoenix, Abell 1835, 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 Mc-Donald 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 Abell 1835, where the [O ii] filaments presented here for the first time are also coincident with molecular gas rising behind the location of X-ray cavities ( Table 1, as well as the continuum predicted by the procedure described in §2.2.1 and Figure 1 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 [O ii] 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 [O ii] 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 §4.2.1.

Star Formation Rates
One of the main goals of our new HST observations was to secure more precise star formation rates (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  Table 4, col. 7), and use the 1st, 2nd, and 3rd quartiles to apply a reddening of E(B − V ) = 0.27 +0.16 −0.07 . Following this intrinsic extinction correction, we also apply a Galactic extinction correction at the redshifted [O ii] wavelengths using the values listed in the NASA/IPAC Extragalactic Database (NED), which are based on Sloan r band data from Schlafly & Finkbeiner (2011). Following these corrections, we converted our corrected [O ii] fluxes to SFRs using the scaling relations found in Kewley et al. (2004): 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 [O ii] 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: where α B = 2.6 × 10 −13 cm 3 s −1 is the coefficient for case B recombination (i.e. only net recombinations, thus excluding transitions directly to the ground state), n e = 300 cm −3 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 log N i = 1.05 log L 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 47 erg s −1 for H1821 (Russell et al. 2010), we find an approximate ionizing photon rate of N i ≈ 10 57 s −1 and a corresponding Strömgren sphere radius of about 0.7 kpc = 0. 16. The amount of quasar contamination to the [O ii] 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 which cannot be included in the SFR calculation. Instead, in lieu of the complexity of modeling the HST PSF, we calculate a radius of 1. 8 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 7 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   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, AGB and HB stars, etc.). 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 (Ṁ cool ; compiled by McDonald et al. 2018) in Figure 3, and can also be found in Table 3 and Table 4 , 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) cool-    Table 3), along with homogenized Hα values from the literature (squares; see Table 4). Curves of constant cooling efficiency ( cool ≡ SF R/Ṁ cool ) are plotted as diagonal black lines, from 1% to 100%. Our best fit relation between SFR andṀ cool obtained from a robust Bayesian analysis is quoted at the top and shown as a solid blue line with shaded 1σ uncertainty band, along with the best fit relations from Fogarty et al. (2015, solid green line) and McDonald et al. (2018, solid red line), and a BCES orthogonal fit (orange line). The steeper-than-unity slope in all these relations suggests that cooling efficiency increases withṀ cool , implying a gradual saturation of AGN feedback, likely tied to an increasing fraction of feedback energy output being radiative as opposed to mechanical. Interestingly, the three starred data points with the highest cooling efficiency are also quasar-hosting systems (Phoenix, H1821, IRAS09104). Points are color coded to correspond to different cooling rate bins, as described in the legend in the top left of this panel. Right: Distributions of cooling efficiency after binning the points from the left panel by cooling rate and using the same color coding. Solid curves overlaid on each histogram are smooth probability density functions calculated non-parametrically via kernel density estimation. It is visually straightforward to see that with higher cooling rate bins, the mean cooling efficiency increases significantly. A slight redshift trend may be due to selection bias, but further analysis is left for a future study.
ing rates and would tend to bringṀ 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. Both SFR andṀ 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Ṁ cool , with much higher cooling efficiencies (defined as the ratio cool = SF R/Ṁ cool ) than the ref-erence 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 wellknown two 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 datapoints 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 datapoints in the scatterplot in Figure 3. Motivated by this trend, we fit the SFR vsṀ cool plot with a total least squares regression using the python software package LinMix (Kelly 2007).
LinMix 4 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 log(SFR) = (1.67 ± 0.17) log(Ṁ cool ) + (−3.25 ± 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Ṁ cool < 100 M yr −1 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 measureṀ 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 log(SFR) = (2.08 ± 0.14) log(Ṁ cool ) + (−4.12 ± 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 ( § 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 quasarmode 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;Sądowski & 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 (Calzadilla et al. in prep).

A Gradual Saturation of AGN Feedback?
The steeper-than-unity relation between SFR and ICM cooling rate presented in § 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 7 K) to warm (10 4 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 non-radiative (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) demonstrate that the accretion rate of SMBHs (Ṁ • ) scaled by the Eddington rate is related to the ICM cooling rate asṀ • /Ṁ Edd ∝Ṁ 1.87 cool , 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 (Ṁ cool 10 3 M yr −1 ), should have 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-Ṁ 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 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Ṁ • /Ṁ Edd → 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;Sądowski & 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 jetinflated 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.Ṁ • 0.1Ṁ Edd , corresponding to a cooling rate of ∼1000 M yr −1 (see Fig. 8 in Mc-Donald et al. 2018). Above this accretion rate is where the quasar-hosting systems in our extreme cooling sample reside (see Fig. 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. Churazov et al. 2001;Li et al. 2017;Reynolds et al. 2002;Yang & Reynolds 2016;Soker et al. 2001;Zhuravleva et al. 2014;Hillel & Soker 2017;Yang et al. 2019;Gaspari 2015). 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 Ṁ • ∼ 0.1Ṁ 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. 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, Mc-Donald 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. 2015Prasad et al. , 2020 which show 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-Ṁ cool plot as well as the decreasing scatter with higher cooling rates shown in Figure 3. For instance, the lack of clusters with bothṀ cool 1000 M yr −1 and SFR 30 M yr −1 , 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. 2015Prasad et al. , 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 (Ṁ cool 1000 M yr −1 ), and thus closer to Eddington accretion rates (Ṁ • 0.1Ṁ 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 ultra-massive black holes up to several 10 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 Myr, 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.

Predicting the Onset of Thermal Instabilities
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 lineof-sight 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. 2010Canning et al. , 2014, with collisional heating being a potential ionization source instead (Ferland et al. 2009). Even so, these Hα and [O ii] maps of Perseus and others still trace where we see warm ∼10 4 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 ten 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 §4.2.2. These extent measurements are all listed in Table 5, and plotted against different Xray 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. 25 th and 75 th ) of the extent distributions, respectively. The larger, dashed, gray y-errorbars on each of these datapoints depicts the minimum and maximum extents across all azimuthal sectors for each cluster, which serves to encode the asymmetry of certain systems like Abell 1795, 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.

R Hα,[Oii] vs Bubbles
One of the mechanisms by which AGN feedback may promote cooling is via radio bubbles inflated by the AGN. This "stimulated feedback" can be achieved by the wake transport of low-entropy warm gas by the buoyantly rising radio bubbles (e.g. McNamara et al. 2016). This uplifted warm gas has an increased infall time, which allows for in-situ production of cold (10 − 100 K) molecular gas. Joint observations of molecular gas with ALMA, and of hot gas with Chandra have revealed the molecular filaments appearing draped around X-ray cavities in a number of systems, with masses and kinematics consistent with uplift by the radio bubbles (e.g. Phoenix: Russell et al. 2017b;Abell 1835: McNamara et al. 2014Abell 1664: Calzadilla et al. 2019. One may expect that from our extreme cooling sample of clusters, we would see similarly intricate networks of filaments preferentially extended along the bubble/jet axis of AGN fueled by the strong cooling flows. However, as we showed in Figure 2, the strongest cooling clusters have a diversity of extended [O ii] nebula morphologies, with some appearing chaotic rather than orderly as in an uplift scenario. For instance, comparing the Phoenix cluster to MACS1931, we see the [O ii] filaments extending mostly along the bubble axes in the former, and perpendicular in the latter. To investigate this point further, we compare our extreme cooling sample to systems in the literature with measured X-ray cavity sizes and cluster-centric distances. Using the sample of Diehl et al. (2008), we collected the cavity size and location measurements (see Table 5) for those clusters that have corresponding Hα extent data from McDonald et al. (2010). In a stimulated feedback scenario, multiphase gas may not extend beyond the maximum radius to which an AGN outburst can uplift low entropy gas. However, in Figure 4  emission (blue points), using maps from McDonald et al. 2010 or this work, respectively) vs different proposed indicators of thermal instability from X-ray observations. In all panels, this average extent was calculated by measuring the maximum extent in ten equally-spaced angular bins, with the y-axis values representing the median over these ten measurements and the (colored) uncertainties representing the interquartile range (i.e. the 25 th and 75 th percentiles). These uncertainties are extended with vertical dotted lines to show the minima and maxima in these measurements. In all panels, a diagonal grey dashed line also represents a one-to-one relationship. Left: R Hα, [Oii] vs extent of X-ray cavities (i.e. radio bubbles) from Diehl et al. (2008), using the distance to the cavity center as well as the leading edge distance for the uncertainties. Middle: R Hα, [Oii] vs the radius where an inflection in the t cool /t ff profile (modeled as two powerlaws) occurs. Right: R Hα, [Oii] vs radius where t cool = 1 Gyr. In middle and right panels, the x-axis values for the blue points come from modeling the ICM profiles in our extreme cooling sample, while values calculated using Hogan et al. (2017) profiles are in red. Uncertainties for each are calculated via bootstrap resampling. In each of these panels, we calculate the Spearman-ρ and Pearson-r correlation coefficients between the two axes. The nearly total asymmetry to one side of the one-to-line in each of these panels demonstrates that all of these predictors of instability establish a volume within which multiphase cooling gas resides.
its leading edge (i.e. cavity center plus cavity radius). The correlation strength is calculated via Pearson-r and Spearman-ρ coefficients, which are suited for assessment of exclusively linear or monotonically correlated data respectively. Both metrics are used to assess the presence and strength of a correlation non-parametrically, with an additional diagnostic of ρ > r potentially indicating a non-linear relationship between two variables. The weak correlation between bubbles and multiphase gas here (ρ = 0.34, p ρ = 0.22; r = 0.36, p r = 0.19) suggests that while bubble uplift is likely a factor in promoting cooling in some systems (e.g. in Perseus and Phoenix), it is not the whole story for most systems.
This weak correlation should perhaps be unsurprising, as gas beneath cavities tends to be turbulent, and will follow the local pressure gradient. Gas filaments that condense out of uplifted low-entropy ICM material should eventually decouple from the bubble wake and fall back down to an altitude where the density contrast is minimized. There will always be such time evolution to the extents of the multiphase gas as well as the cavities that will wash out correlations, which are difficult to account for in observations (e.g. Qiu et al. 2021;Fabian et al. 2021). To complicate matters further, some sys-tems have multiple generations of X-ray cavities (e.g. Perseus: Fabian et al. 2006, Hydra: Wise et al. 2007, NGC5813: Randall et al. 2011, sometimes even perpendicular to each other (e.g. RBS797: Ubertosi et al. 2021), which make it difficult to connect a specific outburst to the extent of cooling seen at longer wavelengths. Conversely, the fact that the star-forming filaments in some systems extend beyond any detected X-ray cavities (e.g. Phoenix, A1795) does not necessarily rule out bubble uplift, but it is impossible to say without deeper data, and even that may not help when older bubbles are projected along the line of sight.
Importantly, however, we can still learn from Figure 4 that filaments on average always lie interior to the radius where we observe bubbles. Thus, the presence and locations of X-ray cavities is still a valuable metric for determining the radius within which the ICM becomes unstable. More measurements of the velocity structure of warm Hα or [O ii] filaments, for instance with integralfield spectroscopy, or of turbulent motion in the ICM in the future (e.g. Hitomi Collaboration et al. 2016; Barret et al. 2016) will help to further extend our understanding of bubble uplift.

R Hα,[Oii] vs 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) 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 chaotic cold accretion 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) argue 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 §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 powerlaws to the t cool /t ff profiles, fixing the slope of just the innermost powerlaw 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 powerlaws cross, and plot these inflection radii (R(t c /t ff ) inflection ) against Hα and [O ii] 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 c /t ff ) inflection . The fact that ρ > r could be an indication of a slightly non-linear relationship between (R(t c /t ff ) 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 . 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) show 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 Abell 2029. We show that there is a very strong correlation between t cool and the average extent of cooling (ρ = 0.88, p ρ = 1.9 × 10 −6 ; r = 0.89, p r = 8.1 × 10 −7 ). The addition of our seven new extent measurements (plotted in blue) are especially helpful in establishing this correlation. Once more, we see that all of the average extents in this 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 datapoints 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 datapoints closest to a one-to-one correspondence with the extent of optical filaments, as shown in Figure 5. This tight relation suggests that, on aver-  Figure 5. Same as in Figure 4, but focusing on the right panel correlation between R Hα, [Oii] and a different t cool threshold. Extended y-errorbars showing the minimum and maximum extents of filaments have been omitted for clarity. We find that measuring the radius at which t cool reaches 0.5 Gyr (rather than 1 Gyr) has an approximately one-to-one correspondence with the extent of filaments.
age, when the ICM reaches a cooling time of 500 Myr or shorter, thermal instabilities will develop, leading to extended filaments of multiphase gas.
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-toone 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 [O ii] 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 maxi-mum 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 7 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.

SUMMARY
In this work, we present new HST observations of the [O ii] emission line nebulae in the strongest starbursting BCGs in the Universe. These new maps allow us to test the limits of AGN feedback in the presence of overwhelming cooling from the ICM. Together with archival Chandra X-ray data, we can link the hot, X-ray emitting phase to this multiphase cooling gas to probe the development of thermal instabilities in the ICM. Our findings can be summarized as follows: This steeper-than-unity relationship means that the cooling of hot gas and the formation of young stars is most efficient in the strongest cool cores.
3. This increasingly efficient conversion of hot (∼ 10 7 K) gas into warm star forming material implies a gradual decrease in the effectiveness of AGN feedback with higher ICM cooling rates. We propose that, as the cooling rate increases, the SMBH accretion rate will approach the Eddington limit, leading to an increasing fraction of the accretion energy released via radiation, rather than via the kinetic mode. The former is less effective at halting large-scale cooling, which would lead to an increase in the global SFR. Under this interpretation, it may not be a coincidence that the most efficiently cooling systems in our sample also host quasars.
4. Using the average extent R of the multiphase gas measured from our [O ii] maps (along with Hα measurements from the literature) as a proxy for where the ICM has become thermally unstable, we compare to features in the ICM to assess how these instabilities develop. We show, 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.
5. We find a strong correlation between R and the cooling radius of the hot ICM. Specifically, we find a one-to-one correlation between the average extent of the multiphase gas and the radius at which the ICM cooling time reaches 0.5 Gyr, which may be indicative of a universal condensation timescale in cluster cores.
The new data presented here represent the sharpest view yet of the massive star forming regions in the strongest starbursting BCGs in the Universe. The unique environments provided by these systems allow us to test the limits of AGN feedback in the presence of overwhelming cooling. These systems could possibly mimic the environments of higher redshift cluster cores, as there is evidence that they are more likely to harbor central quasars as well as starbursts (e.g. Hlavacek-Larrondo et al. 2013;Somboonpanyakul et al. 2022), though possibly fueled somewhat differently (e.g. Mc-Donald et al. 2016). Thus, this extreme cooling sample offers us a low-redshift window into higher-redshift phe-nomena, making these ideal candidates for future followup with observatories like the James Webb Space Telescope.

ACKNOWLEDGMENTS
This work is based on observations with the NASA/ESA Hubble Space Telescope obtained at the Space Telescope Science Institute, which is operated by the Associations of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. These observations are associated with programs 15315, 15661, and 16001. MC acknowledges support from the NASA Headquarters under the Future Investigators in NASA Earth and Space Science and Technology (FINESST) award 20-Astro20-0037. MM and MC acknowledge financial support from programs HST-GO15315, HST-GO15661, and HST-GO-16001, which was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Associations of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. MM and MC acknowledge additional financial support for this work, provided by the National Aeronautics and Space Administration through Chandra Award Number GO0-21114A issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. MG acknowledges partial support by NASA Chandra GO9-20114X and HST GO-15890.020/023-A, and the Black-HoleWeather program.
This work has made extensive use of the SAO/NASA Astrophysics Data System (ADS) and the arXiv preprint server.