The Relationship Between Brightest Cluster Galaxy Star Formation and the Intracluster Medium in CLASH

A detailed summary of the science-level data products for CLASH may be found in Postman et al. (2012). For the 11 X-ray -elected clusters with evidence of BCG star formation activity, we used the CLASH photometric data set covering 16 filters spanning from ∼2000 to 17000 ${\rm{\mathring{\rm A} }}$ in the observer frame. We used drizzled mosaics with a 65 milliarcsecond pixel scale, the same image data used in Fogarty et al. (2015). These data are publicly available via MAST HLSP.⁷ Drizzled mosaics were constructed using the MosaicDrizzle pipeline (Koekemoer et al. 2011). In keeping with our previous work, we calculated a single Milky Way foreground reddening correction for each BCG in each filter using the Schlegel et al. (1998) dust maps. Our drizzled images are background corrected using an iterative 3σ clipping technique. Owing to additional uncertainty in the flat-fielding of WFC3/UVIS photometry over large spatial scales, discussed in Fogarty et al. (2015), the median flux measured in an annulus around each BCG was subtracted from each of the UV filters as well.

Spitzer/IRAC 3.6 μm and 4.5 μm mosaicked observations and catalogs are available for all of the BCGs studied in this paper, and 5.7 μm and 7.9 μm observations are available for Abell 383, MACS 1423.8+2404, and RXJ 1347.5−1145. Fluxes for Spitzer IR sources in the CLASH fields were taken from the publicly available CLASH/Spitzer catalog.⁸ We describe aperture selection and our method for correcting for crowded fields in the Spitzer photometry in Section 3.1.

Photometry for each channel was measured on mosaic images generated using the MOPEX software package. The default MOPEX settings for the catalog use the 3.6 μm channel, and use the fiducial image frame for this channel to generate the mosaics for longer wavelength IRAC channels for each CLASH cluster. The Spitzer mosaic images have a pixel scale of 0.6 arcsec. Flux values were obtained with Source Extractor in double-image mode using MOPEX weights (Bertin & Arnouts 1996). The Source Extractor-generated catalogs consist of aperture photometry for sources in fixed-diameter apertures ranging from 2 to 40 pixels. Full details about the parameters used with MOPEX and Source Extractor to generate the CLASH catalog are given in the online documentation.⁹

We used archival data from Herschel/PACS (100, 160 μm) and Herschel/SPIRE (250, 350, and 500 μm) to extend BCG photometry into the far-IR. Table 1 details the observations used and their exposure times. We obtained level 2 archival data using the Herschel Science Archive reduction pipeline for SPIRE observations and the HSA MADMAP reduction pipeline for PACS observations. Observations were co-added and photometric parameters were measured using the HIPE software package. We describe Herschel aperture selection and background subtraction in Section 3.1.

Chandra data exist for all CLASH clusters, with exposure times for individual observations ranging from 19.3 ks (for Abell 383) to 115.1 ks (for MACS 1423.8+2404). These data were assembled in Donahue et al. (2014) to construct the parameter profiles that we use to analyze the thermodynamics of the ICM in the environs of each BCG. We also examine the archival data for evidence of X-ray-loud AGNs, in order to determine whether it is important to consider AGN emission models in our SED fits. We found evidence for only one X-ray-loud AGN, in the cluster MACS 1931.8−2635. The AGN classification is based on the presence of a Chandra point source, with a spectrum showing evidence for hard X-ray (>5 keV) emission over that expected from hot gas with ${kT}\sim 5\mbox{--}7\,\mathrm{keV}$.

We constructed multi-instrument UV through IR SEDs for the UV-bright filamentary features in each UV-luminous CLASH BCG. Since we are interested in measuring the star formation properties of the BCGs, for the Hubble UV through near-IR photometry, we extracted flux from within the apertures described in Fogarty et al. (2015; apertures for each BCG are shown in Figure 3 in that paper). These apertures were selected to encompass the region in each BCG exhibiting a UV luminosity of $7.14\times {10}^{24}$ erg s⁻¹ Hz⁻¹ pix⁻², corresponding to a star formation surface density of ${{\rm{\Sigma }}}_{\mathrm{SFR}}\geqslant 0.001\,{M}_{\odot }$ yr⁻¹ pix⁻², after accounting for dust reddening. Measuring fluxes in these apertures maximizes the contribution made by the recently formed stellar population to the SED, while it minimizes the contribution made by the passive stellar population in the bulk of the BCGs. Our procedure for measuring fluxes minimizes the risk of underestimating dust attenuation in the star-forming regions of the BCGs since we are not averaging dust attenuation over the dusty star-forming and relatively dust-free quiescent parts of the galaxy. We do not match the Hubble point-spread functions (PSF) since the photometric aperture sizes we use are much larger than the sizes of the PSFs for the HST passbands.

Estimating mid-IR fluxes using Spitzer was more complicated, since star-forming features are not well-resolved spatially in Spitzer photometry. We measured IR fluxes from Spitzer in apertures in the CLASH/Spitzer catalogs that were selected to encompass the BCG while minimizing inclusion of satellite galaxies. Both the old stellar population and dust emission (primarily in the form of PAH features) contribute significantly to the flux in these filters. Since the angular resolution of Spitzer (>1.45 arcsec) is insufficient to resolve the star-forming features we wish to study, we needed to subtract the contribution made by old stellar light from outside our Hubble apertures to the Spitzer IR fluxes. We accomplished this by first measuring the HST/WFC3 F160W flux in the region covered by the Spitzer apertures but not the Hubble apertures. We assumed that this flux is due entirely to old stellar light, and estimated the (IRAC − F160W) color by modeling old stellar populations experiencing dust attenuation $\leqslant 0.5\,{A}_{{\rm{V}}}$ using the Bayesian SED-fitting algorithm iSEDfit (see Section 3.2 for details) for each BCG. We used these colors to scale the F160W fluxes corresponding to the old stellar population outside the Hubble apertures to Spitzer IRAC fluxes and subtracted these scaled fluxes from our Spitzer photometry. The resulting Spitzer IRAC photometry corresponds to the fluxes in the Hubble apertures. An example pair of Spitzer and Hubble apertures is shown in Figure 1.

Zoom In Zoom Out Reset image size

We note that ${{\rm{H}}}_{2}$ vibrational modes and other near-IR emission lines have been detected between 5 and 25 μm in star-forming BCGs (Donahue et al. 2011). However, contamination from emission lines similar to those observed in Donahue et al. (2011) would only contribute a few percent to the mid-IR fluxes derived for the stellar and dust emission in star-forming regions, while the uncertainties on these fluxes after subtracting the excess old stellar light component is $\gtrsim 20 \%$. Therefore, we do not attempt to estimate the contribution of these lines to the Spitzer mid-IR photometry.

We extracted Herschel PACS and SPIRE photometry from archival data using the HIPE software package. Since Herschel bands are dominated by dust re-emission of UV flux, we assumed that all of the Herschel flux in each BCG comes from the star-forming features. We measured photometry in circular apertures large enough to cover the PSF of each Herschel filter. Aperture radii were chosen to be 8'', 12'', 18'', 24'', and 36'' for the 100 μm, 160 μm, 250 μm, 350 μm, and 500 μm filters, respectively. To obtain fluxes, we took the mean of pairs of cross scans, while we estimated uncertainties using the difference between scans. For each Herschel band, we found typical uncertainties on the order of $\lesssim 2 \%$ for the PACS filters and $\sim 10 \%$ for the SPIRE filters. Since the BCGs occupy crowded fields, we measured crowded source backgrounds in annuli centered on the BCG with inner radii of 16'', 24'', 36'', 48'', and 72'' for the respective PACS and SPIRE filters and outer radii of 180''.

3.1.1. Photometric Errors

We fit SEDs for the CLASH BCGs using the Bayesian fitting code iSEDfit (Moustakas et al. 2013). In order to take full advantage of the CLASH UV–FIR observations, we incorporated modifications to iSEDfit to account for emission from dust. The iSEDfit package fits SEDs by first generating a grid of parameters that describe models of dust-attenuated stellar emission by a composite stellar population. Grids are produced by randomly sampling a bounded section of parameter space. The sampling is weighted with either a uniform prior or a logarithmic scale prior. With these parameters, iSEDfit constructs a grid of model spectra given a choice of synthetic stellar population (SSP), initial mass function (IMF), and dust-attenuation model. The package then computes the synthetic photometry by convolving these model spectra with the filter responses of the CLASH SEDs and uses this grid of synthetic SEDs to sample the posterior probability distribution that the model stellar population and dust parameters produce the observed SED.

The iSEDfit package uses a Bayesian Monte Carlo approach to estimate the probability distributions for model parameters. For each BCG, we constructed a model grid consisting of 20,000 models. The likelihood and mass scaling ${ \mathcal A }$ of each model given the observed SED was calculated by minimizing ${\chi }^{2}$ for each model,

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{({F}_{i}-{ \mathcal A }{M}_{i})}^{2}}{{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 1 }$$

where F_i and ${\sigma }_{i}$ are the SED fluxes and uncertainties, and M_i are the model fluxes. The minimum ${\chi }^{2}$ and ${ \mathcal A }$ are found by solving for $\partial {\chi }^{2}/\partial { \mathcal A }=0$. We obtained a posterior probability distribution for the models using a weighted random sampling of the model grid, where model weights were determined by the likelihoods. Posterior probability distributions for individual physical parameters were obtained by taking the distribution of parameters of the sampled models. For the SFR, the posterior consists of the distribution of the instantaneous normalized SFR determined by the parameterized star formation history, multiplied by ${ \mathcal A }$. For a detailed discussion of iSEDfit, see Moustakas et al. (2013).

We model the star formation history of each BCG with an exponentially decaying curve (the early-type population) and a superimposed exponentially decaying burst at recent times. The initial exponentially decaying curve is parameterized by the age of the BCG t and the decay rate of the curve τ, while the exponentially decaying burst is parameterized by the duration of the starburst ${\rm{\Delta }}{t}_{b}$, the decay rate of the burst, and the mass of the burst relative to the old stellar population (see Figure 2). The bounds on parameter space for the entire model and the assumptions we made for the BCG stellar populations are given in Table 2.

Zoom In Zoom Out Reset image size

Table 2.  SED-fitting Parameters

Stellar Population Model
Synthetic Stellar Population	Bruzual & Charlot (2003)
Initial Mass Function	Salpeter (1955)
Attenuation Law	Calzetti et al. (2000)
Dust Emission	Draine & Li (2007)

	Minimum	Maximum	Sampling
Model Parameter Space Constraints	Value	Value	Interval
Old Stellar Population
Age, t	6 Gyr	zAgea	Linear
Decay Rate, τ	$0.05t$	$0.2t$	Linear
Metallicity	$3\times {10}^{-2}\,{Z}_{\odot }$	$1.5\,{Z}_{\odot }$	Linear
Burst Population
Burst Duration, ${\rm{\Delta }}{t}_{b}$	10−2 Gyr	5.0 Gyr	Logarithmicb
Burst Decay Percentage	0.01	0.99	Linear
Relative Burst Mass, fburstc	0.0016	6.4	Logarithmic
Dust Parameters
Attenuation ${A}_{V}$	0	2	Linear
PAH Abundance Index qPAH	0.10	4.58	Lineard
${\gamma }^{e}$	0.0	1.0	Linear
Umine	0.10	25.0	Logarithmic
Umaxe	103	107	Logarithmic

Notes.

^aZ_Age is the age of the universe at the BCG redshift Z. ^bBurst parameters were sampled logarithmically, since their qualitative effect on the model SED of the galaxy occurs on order-of-magnitude scales. The exception to this is the burst decay percentage, which is one minus the amplitude of current star formation activity relative to the amplitude of the burst ${\rm{\Delta }}{t}_{b}$ yr ago. ^cMass of stars created by the starburst at t relative to the mass of stars created by the exponentially decaying old stellar population at t. The burst mass percentage is calculated using fburst/(1+fburst). ^dThe Draine & Li (2007) model parameters sampling intervals were chosen based on the model parameter distributions of the template spectra. ^eDraine & Li (2007) treat dust in a galaxy as consisting of two components. The first component consists of a fraction γ of the dust that is exposed to a power-law distribution of starlight intensity, ranging from U_min to U_max, while the second component consists of the remainder of the dust and is only exposed to a starlight intensity ${U}_{\min }.\ U$ is defined to be the intensity of starlight relative to the local radiation field, and U_min and U_max are bounds on the distribution of U.

Download table as:  ASCIITypeset image

We modified the SED-fitting routine to incorporate dust emission models from Draine & Li (2007) into the parameter grid and synthetic spectra used in iSEDfit. The dust parameter space was sampled using the bounds and priors given in Table 2. For each synthetic spectrum, the total luminosity of the dust spectrum was normalized by the difference between the unattenuated and attenuated stellar spectrum. Figure 3 shows an example of the synthetic attenuated stellar plus dust spectrum, and the Appendix shows best-fit synthetic spectrum for each BCG.

Zoom In Zoom Out Reset image size

With our modified version of iSEDfit, we were able to test multiple stellar population and dust-attenuation models in order to examine the dependence of our results on the choice of SSP, IMF, and attenuation law. We adopt a Bruzual & Charlot (2003) SSP, the Salpeter (1955) IMF defined over the interval 0.1–100 ${M}_{\odot }$, and the Calzetti et al. (2000) attenuation law. We chose to adopt a Salpeter (1955) IMF in order to produce SFRs consistent with the SFRs estimated in our previous paper using the Kennicutt (1998) relationship. In order to examine the effect of an ${A}_{{\rm{V}}}$-dependent attenuation curve, we also performed fits using a modified Calzetti attenuation law, where we multiplied the Calzetti curve with an attenuation-dependent slope that matches the attenuation dependence of the curves published in Witt & Gordon (2000). To test if our results have any dependence on our choice to adopt a Calzetti attenuation law, we also ran iSEDfit assuming clumpy SMC-like dust in a shell geometry. We found our results to be largely consistent with results using the regular Calzetti curve, as well as the attenuation-dependent version of the Calzetti curve. This is because, in the typical CLASH BCG, the attenuation is ${A}_{{\rm{V}}}\lesssim 1.0$, where the overall attenuation dependence on the shape of the curve only has a modest effect. As a further test of the model dependence of our results, we performed SED fits assuming a Chabrier (2003) IMF. Altering the IMF shifts the SFRs downwards by 0.25 dex, but otherwise does not significantly alter the stellar and dust parameters we seek to measure (${\rm{\Delta }}{t}_{b}$, ${M}_{d},{A}_{{\rm{V}}}$).

We calculated the radial profiles of the cooling time, defined as the ratio of the thermal energy density to the rate of radiative energy density loss for an optically thin plasma undergoing Bremsstrahlung emission, for each cluster using

$$\begin{eqnarray}{t}_{\mathrm{cool}}(r) & = & \displaystyle \frac{3}{2}\displaystyle \frac{n(r){k}_{{\rm{B}}}T(r)}{{n}_{{\rm{e}}}(r){n}_{{\rm{H}}}(r){\rm{\Lambda }}[T(r),Z(r)]}\\ & = & \displaystyle \frac{6.9}{2}\displaystyle \frac{{k}_{{\rm{B}}}T(r)}{{n}_{{\rm{e}}}(r){\rm{\Lambda }}[T(r),Z(r)]},\end{eqnarray} \tag{ 2 }$$

where k_B is the Boltzmann constant, $n(r)$ is the total number density profile of the plasma, ${n}_{{\rm{H}}}(r)$ is the H number density, ${n}_{{\rm{e}}}(r)$ is the electron number density, $T(r)$ is the temperature profile, $Z(r)$ the metallicity profile, and ${\rm{\Lambda }}(T,Z)$ is the cooling function. We obtained values of the cooling function for specific temperatures and metallicities by interpolating the cooling function values in Sutherland & Dopita (1993), and assumed $n\approx 2.3{n}_{{\rm{H}}}$ (Cavagnolo et al. 2009). The ICM density, temperature, and metallicity profiles used in this study are available in Donahue et al. (2014). We used the non-parametric Joint Analysis of Cluster Observations (JACO) profiles (see Mahdavi et al. 2007 and Mahdavi et al. 2013 for a description of the JACO algorithm), which are reported in concentric shells spaced so that each annulus contains at least 1500 X-ray counts.

We measured t_cool at specific radii by interpolating on ${n}_{{\rm{e}}}(r),T(r)$, and $Z(r)$ at the desired radius and solving Equation (2). In order to determine the uncertainty on t_cool, we produced an ensemble of 1000 ${n}_{{\rm{e}}},T$, and Z profiles, where the values in each profile were obtained by drawing from normal distributions defined by the observed values and uncertainties in each profile. By calculating the distribution of t_cool given the interpolated values of ${n}_{{\rm{e}}},T$, and Z for the ensemble of profiles, we sampled the probability distribution of t_cool and thus estimated the uncertainties. We constrained the ensemble of metallicity profiles to only allow values of Z in the range $0.15\mbox{--}1.5{Z}_{\odot }$, since ICM metallicities are typically $0.3{Z}_{\odot }$ and tend to be $\gtrsim 0.6\mbox{--}0.8{Z}_{\odot }$ in the centers of cool-core clusters (De Grandi et al. 2004).

We calculated freefall times by assuming the cluster masses obey an NFW profile (Navarro et al. 1997). NFW mass concentration parameters and values of ${M}_{200}$ were obtained from Merten et al. (2015). We computed the enclosed mass for each cluster as a function of radius,

$$\begin{eqnarray}&&{M}_{\mathrm{enc}}(r)=4\pi {\rho }_{0}{r}_{s}^{3}\left[\mathrm{ln}\left(\displaystyle \frac{{r}_{s}+r}{{r}_{s}}\right)-\displaystyle \frac{r}{{r}_{s}+r}\right],\end{eqnarray} \tag{ 3 }$$

where ${\rho }_{0}$ and r_s are NFW scale factors determined by the mass and concentration parameter of the galaxy cluster, as well as the critical density at the cluster redshift. We then calculated freefall time profiles,

$$\begin{eqnarray}&&{t}_{\mathrm{ff}}(r)=\sqrt{\displaystyle \frac{2{r}^{3}}{{{GM}}_{\mathrm{enc}}(r)}}.\end{eqnarray} \tag{ 4 }$$

We estimated the contribution of the BCG stellar mass to the freefall time using the Cooke et al. (2016) stellar mass estimates of CLASH BCGs. We did not estimate BCG stellar masses with our SED fits since our photometry does not cover the entire extent of the BCG. Stellar mass density profiles were estimated from the stellar mass of each BCG assuming a Hernquist (1990) profile and the Shen et al. (2003) mass−size relationship (Ruszkowski & Springel 2009; Laporte et al. 2013). The inclusion of stellar mass in our estimates of ${M}_{\mathrm{enc}}(r)$ does not significantly affect the freefall time profiles at radii $\gtrsim 25\,\mathrm{kpc}$ (or similarly $\gtrsim 0.025\,{R}_{500}$), which is not surprising given that the BCG stellar mass dominates the cluster density profile only within the central ∼10 kpc of a cluster (Newman et al. 2013; Caminha et al. 2017; Monna et al. 2017). Still, we include the BCG stellar mass in our estimation of the freefall time profiles and find that the BCG stellar masses alter the freefall times by $\lesssim 4 \%$ at ∼10 kpc and $\lesssim 1 \%$ at ∼20 kpc.

Comparing the characteristics of BCG starbursts to ICM thermodynamic parameters requires choosing a radius in the ICM profile at which to measure these parameters. We measure ICM parameters at $0.025\,{R}_{500}$, which is the smallest fraction of R₅₀₀ that does not require extrapolation of the X-ray profiles for the 11 CLASH clusters studied. We chose this radius since we expect feedback effects to be strongest near the centers of cool cores and negligible outside the core. At radii outside the cool core, the thermodynamical state of the ICM is only weakly tied to feedback and cooling, and therefore should be weakly related to the BCG (Cowie et al. 1983; Hu et al. 1985).

We also base our expectation of the radial dependence between the properties of the ICM and feedback on the BCG simulations in Li et al. (2015). At large radii ($\gtrsim 100$ kpc), the variation in ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ with time in the simulation (and with the level of star formation and AGN-driven feedback) is weak, so a relationship between SFR and cluster dynamical state would be difficult to detect (see their Figure 9). On the other hand, ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ varies wildly at small radii (≲10 kpc), so in this case, too, relationships involving ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ and other parameters related to cooling and feedback in the cluster would be difficult to interpret (although the variation in simulations may be due to the idealized AGN entering periods of complete quiescence as the fuel supply reaches zero). However, since $0.025\,{R}_{500}$ is greater than 10 kpc for all of the CLASH clusters, even if this effect is physically realistic we would not expect to observe it.

For BCGs with detectable star formation activity, we find a tight correlation between the BCG SFR and the ratio ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ measured at $0.025\,{R}_{500}$. The trend and data are shown in Figure 4. The best-fit line between these two parameters is

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\displaystyle \frac{{t}_{\mathrm{cool}}}{{t}_{\mathrm{ff}}}=(1.6\pm 0.4)-(0.15\pm 0.03){\mathrm{log}}_{10}\displaystyle \frac{{\rm{SFR}}}{{M}_{\odot }\,{{\rm{yr}}}^{-1}},\end{eqnarray} \tag{ 6 }$$

and is shown as the black curve in Figure 4. The data are fit by this trend with no detectable intrinsic scatter (${\sigma }_{i}\lt 0.15$ dex at 3σ). We estimated the best-fit lines and intrinsic scatters using the least-squares method in Hogg et al. (2010) for fitting data with uncertainties in two dimensions and intrinsic scatter. The data have a Spearman correlation coefficient ${R}_{S}=-0.98$ and a Pearson correlation coefficient $R=-0.95$.

Zoom In Zoom Out Reset image size

We estimated both the strength of the relationship we observe and the strength of the claim that the relationship has little intrinsic scatter. We first investigated the possibility that there is no relationship between SFR and ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ by calculating the marginal probability distribution of the slope using Equation (35) of Hogg et al. (2010), and find that the slope of the relationship is less than 0 at 99.993% confidence.

We then examined the possibility that a relationship obeying the trend in Equation (6) but with a large intrinsic scatter produced the results we observe. We created an ensemble with 10⁴ synthetic data sets, consisting of 11 SFR values randomly sampled in logarithmic units from 0.1 ${M}_{\odot }$ yr⁻¹ to 1000 ${M}_{\odot }$ yr⁻¹, and calculated the corresponding predicted values of ${\mathrm{log}}_{10}{t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$. Assuming a given intrinsic scatter, ${\sigma }_{i}$, we generated an offset from the trendline for each of the 11 points by sampling a normal distribution of standard deviation ${\sigma }_{i}$ to obtain the amplitude of the offset and determined a direction of the offset by sampling a uniform distribution between 0 and $2\pi$. We added an additional offset for each point in the x-direction and in the y-direction by sampling normal distributions with standard deviations equal to the mean uncertainty of the observed SFRs and the observed ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ values, respectively. The probability that data more tightly correlated than our data set ($\parallel R\parallel =0.95$) would be produced given ${\sigma }_{i}$ was calculated for the ensemble of data sets. Since we used log-normally distributed variables to generate our synthetic data, we used the Pearson coefficient R to measure the strength of the correlation. We repeated this procedure for ${\sigma }_{i}$ ranging from 0.0 to 0.6, and found that for ${\sigma }_{i}\gtrsim 0.22$ dex, the probability of drawing a data set with a trend tighter than the one we observe is $\lt 0.3 \%$. The results of the above correlation test are presented in Figure 5.

Zoom In Zoom Out Reset image size

The ratio ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ is thought to be a proxy for the thermal instability of ICM gas—when the gas can cool quickly relative to the time it takes for it to infall, it can more readily become thermally unstable and collapse (Gaspari et al. 2012; Li & Bryan 2014). Assuming ICM instability increases with decreasing ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$, our results suggest that BCG SFRs will increase with the thermal instability of the surrounding ICM.

The relationship between the star formation duration ${\rm{\Delta }}{t}_{b}$ and cooling time at $0.025\,{R}_{500},{t}_{{\rm{cool}}}$ is shown in Figure 6. The relationship between these two parameters is difficult to quantify owing to the large uncertainties on ${\rm{\Delta }}{t}_{b}$. In order to estimate the strength of this correlation, we calculated the Spearman coefficients both including and excluding MACS 1931.8−2653. In performing these computations, we also need to address the fact that some of the ${\rm{\Delta }}{t}_{b}$ values are only lower limits. If we use the lower limit values to compute the ${\rm{\Delta }}{t}_{b}$ ranking, then ${\rm{\Delta }}{t}_{b}$ and t_cool are strongly correlated, with a Spearman coefficient of 0.86 including MACS 1931.8−2653, and a coefficient of 0.81 excluding it. If we assume that the lower limits in ${\rm{\Delta }}{t}_{b}$ are tied for the highest rank, then the Spearman coefficient including MACS 1931.8−2653 is 0.57, implying a weak correlation ($p\lt 0.065$), while if MACS 1931.8−2653 is excluded, the data are not significantly correlated (${R}_{S}=0.427,p\lt 0.22$). In the limit of the worst possible ${\rm{\Delta }}{t}_{b}$ ranking for the lower limit data (e.g., the data point with the highest ranked t_cool among the "lower limit" points has the lowest relative ${\rm{\Delta }}{t}_{b}$ rank), the two data sets are not correlated (${R}_{S}=0.433,p\lt 0.18$ with MACS 1931.8−2653, ${R}_{S}=0.24,p\lt 0.50$ without). Therefore, cooling times and burst durations are likely weakly positively correlated (although this observation may be driven by MACS 1931.8−2653), and ${\rm{\Delta }}{t}_{b}$ and ${t}_{{\rm{cool}}}$ become comparable in amplitude when ${\rm{\Delta }}{t}_{b}$ reaches gigayear timescales.

Zoom In Zoom Out Reset image size

All of the BCGs either lie to the left of the line where ${\rm{\Delta }}{t}_{b}={t}_{{\rm{cool}}}$ or are consistent with lying to the left of this line. This is the region where starbursts occur on timescales shorter than the timescale for the ICM to cool. In clusters where ${\rm{\Delta }}{t}_{b}$ and t_cool fall to the left of this line, star formation has been ongoing for a shorter duration than the time it would take for the ICM at $0.025\,{R}_{500}$ to radiatively cool. Therefore, assuming t_cool profiles typically increase monotonically with radius, if the cold gas reservoir fueling star formation is the result of the ICM radiatively cooling, then the body of low-cooling-time gas inside $0.025\,{R}_{500}$ that was present at the onset of star formation will be depleted. Alternatively, in clusters where ${\rm{\Delta }}{t}_{b}$ and t_cool fall to the right of this line, either higher-cooling-time gas will have radiatively cooled and replenished the gas inside $0.025\,{R}_{500}$ in order to continue forming the cold gas fueling star formation, or radiative cooling has been arrested and t_cool is locally static.

We find that both the dust mass and burst duration of CLASH BCGs are correlated with their SFRs. The relationship between ${\rm{\Delta }}{t}_{b}$ and SFR is shown in Figure 7. Large SFRs are consistent with star formation episodes that have recently begun, and as the bursts persist to ∼Gyr timescales, the SFRs diminish by several orders of magnitude.

Zoom In Zoom Out Reset image size

Figure 8 shows the relationship between M_d and SFR. In order to be consistent with studies of SFR and ${M}_{d}$ conducted by da Cunha et al. (2010) and Hjorth et al. (2014), when comparing these two quantities, we analyze SFRs and dust masses obtained with a Chabrier (2003) IMF (see Table 5). We continue to use Salpeter (1955) when discussing results in the rest of our paper.

Zoom In Zoom Out Reset image size

Dust masses and SFRs are correlated, with a best-fit trendline of

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\displaystyle \frac{{M}_{d}}{{M}_{\odot }}=({7.1}_{-0.3}^{+0.4})+({0.97}_{0.24}^{+0.34}){\mathrm{log}}_{10}\displaystyle \frac{{\rm{SFR}}}{{M}_{\odot }\,{{\rm{yr}}}^{-1}},\end{eqnarray} \tag{ 7 }$$

with an intrinsic scatter of $\lt 0.83$ dex ($3\sigma$ limit). The analysis of da Cunha et al. (2010) derived dust masses and SFRs for field galaxies using a UV–IR SED-fitting technique that was similar to ours. We overlay their best-fit trendline, which has a best-fit slope of 1.11 ± 0.01 and intercept of 7.1 ± 0.01 in Figure 8. In order to constrain the behavior of the trend at the large SFR end, we include the dust mass and SFR of the Phoenix Cluster BCG (e.g., McDonald et al. 2013). We adopt the Mittal et al. (2017) estimate of 454–494 ${M}_{\odot }$ yr⁻¹ for the Phoenix SFR and fit the far-IR SED of Phoenix to estimate M_d. This extended data set shows a flattening slope at the high-SFR end of the relationship. When we overlay a trend drawn from Hjorth et al. (2014) that takes into account the evolution of star formation and dust for starbursting galaxies with large ($\gtrsim 1000\,{M}_{\odot }$ yr⁻¹) SFRs, we find excellent correspondence with our data across the range of SFRs studied. Although our confidence that we observe a change in the SFR–${M}_{d}$ trend at the high-SFR end is limited by the size of our sample, our results are fully consistent with star-forming BCGs producing dust reservoirs like starbursts in the field.

We wish to understand the possible impact of AGN emission on our SED-fitting results. Therefore, we analyzed the impact of AGN emission on the SED fit to MACS 1931.8−2653, since this BCG shows the strongest evidence for AGN emission. If the AGN component has little impact on the SED fit to MACS 1931.8−2653, we do not believe there will be substantial AGN contamination of our fits to the other starbursting BCGs in CLASH, which host weaker or quiescent AGNs. The X-ray point source in this BCG is well-fit by an AGN power spectrum with a 2–10 keV luminosity of ${5.32}_{-0.37}^{+0.40}\times {10}^{43}$ ergs s⁻¹ cm⁻² (Hlavacek-Larrondo et al. 2013). Moreover, the IR template-fitting results of Santos et al. (2016) imply that a substantial AGN contribution may exist in the case of MACS 1931.8−2653, so it is important to investigate the effect of adding an AGN component to the model used to fit the UV–IR SED of this X-ray-loud BCG.

As Santos et al. (2016) postulate a potentially large AGN contribution to the UV–IR SED of MACS 1931.8−2653, we ran a separate SED-fitting analysis for this BCG wherein we include the IR AGN model of Siebenmorgen et al. (2015). We allowed iSEDfit to sample the full range of parameters that describe the Siebenmorgen et al. (2015) AGN emission model library, which for our purposes are nuisance parameters. The contribution of the AGN to the total UV–IR luminosity of the galaxy was allowed to range between $0.001\times$ and $10\times$ the stellar contribution.

We find that the effect of AGN emission in MACS 1931.8−2653 is marginal. The best-fit ${\chi }^{2}$ degrades slightly, to 1.09, most likely owing to the addition of AGN model nuisance parameters. After incorporating AGN emission, we find the ${\mathrm{log}}_{10}$ $\mathrm{SFR}=2.36\pm 0.19\,{M}_{\odot }$ yr⁻¹, ${\mathrm{log}}_{10}$ ${\rm{\Delta }}{t}_{b}=-{0.69}_{-0.80}^{+0.46}$ Gyr, and that ${\mathrm{log}}_{10}$ ${M}_{d}={8.86}_{-0.38}^{+0.34}\,{M}_{\odot }$. Compared to our AGN-free model, we find that the SFR changes by $\sim 0.3\sigma ,{\rm{\Delta }}{t}_{b}$ by $\sim 0.5\sigma$, and ${M}_{d}$ by $\sim 0.05\sigma$. This is not surprising, since the AGN contribution to the UV–FIR luminosity relative to the stellar contribution is ${\mathrm{log}}_{10}$ ${f}_{\mathrm{AGN}}=-1.65\pm 0.87$.

One possible explanation for the discrepancy between these findings and the results published in Santos et al. (2016) for MACS 1931.8−2653 is our decision not to incorporate WISE photometry in our SED fitting. WISE W3 and W4 filters cover the region of the MIR spectrum (∼10 μm) sensitive to the contribution of AGN flux. Since for our purposes the AGN flux is a contaminant, it does not make sense to include these filters in our SED. The difference may also be explained, in part, by the fitting technique—the theoretical modeling we employ to fit the data allows us greater flexibility to fit the data than the empirical templates used in Santos et al. (2016).

Models of AGN-driven condensation and precipitation that involve ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ as a proxy for thermal instability in the ICM provide a natural foundation for interpreting our observed correlation between ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ and the SFR. If ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ is related to the rate of molecular gas production around an active BCG, it should be correlated with the SFR. Specifically, recent simulations show that ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ determines the critical overdensity at which ICM density perturbations can condense (Singh & Sharma 2015). We hypothesize this critical overdensity, in turn, determines the mass deposition efficiency ${\epsilon }_{\mathrm{MD}}$, defined to be the mass fraction in a region of the ICM that cools into molecular gas. Assuming that density perturbations of the ICM follow a log-normal distribution in the core of the cluster, ${\epsilon }_{\mathrm{MD}}$ at a radius r ought to be

$$\begin{eqnarray}{\epsilon }_{\mathrm{MD}} & = & {\displaystyle \int }_{{\delta }_{c}\left(\tfrac{{t}_{\mathrm{cool}}}{{t}_{\mathrm{ff}}}\right)}^{\infty }\displaystyle \frac{1}{\sqrt{2\pi }\sigma \rho }{e}^{-\tfrac{{(\mathrm{ln}\rho )}^{2}}{2{\sigma }^{2}}}d\rho \\ & = & \displaystyle \frac{1}{2}{Erfc}\left(\displaystyle \frac{\mathrm{ln}{\delta }_{c}\left(\tfrac{{t}_{\mathrm{cool}}}{{t}_{\mathrm{ff}}}\right)}{\sqrt{2}\sigma }\right),\end{eqnarray} \tag{ 8 }$$

where ${\delta }_{c}\left(\tfrac{{t}_{\mathrm{cool}}}{{t}_{\mathrm{ff}}}\right)$ is the critical overdensity for ICM condensation as a function of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ and σ is the width of the ICM density perturbation distribution. The observed relation between SFR and ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ would therefore imply that SFR scales with ${\epsilon }_{\mathrm{MD}}$.

We can use the ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$–SFR relationship to infer properties of the relationship between ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ and ${\epsilon }_{\mathrm{MD}}$, and therefore constrain models of feedback-regulated cooling. As a first-order approximation, we assume that a condition close to equilibrium exists between cooling and star formation for most of the duration of an episode of feedback-regulated cooling, so that SFR $\sim \,{\dot{M}}_{g,\mathrm{real}}$, where ${\dot{M}}_{g,\mathrm{real}}$ is the actual ICM cooling rate. By measuring ${\dot{M}}_{g}\equiv {M}_{g}/{t}_{\mathrm{cool}}$ within $0.025\,{R}_{500}$ with the X-ray parameter profiles used in this paper, we find that

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\epsilon }_{\mathrm{MD}}\sim {\mathrm{log}}_{10}{\rm{SFR}}-{\mathrm{log}}_{10}{\dot{M}}_{g}=\displaystyle \frac{{t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}-10}{-{17.1}_{-2.3}^{+1.9}},\end{eqnarray} \tag{ 9 }$$

where for the purposes of this simple model we have set ${\mathrm{log}}_{10}{\epsilon }_{\mathrm{MD}}=0.0$ when ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}=10.0$, which is a "critical" value for mechanical feedback-triggered condensation (Voit & Donahue 2015). In a relatively uniform sample of clusters like CLASH, where core gas masses occupy a narrow range (∼0.4–1.4 $\times \,{10}^{13}\,{M}_{\odot }$), ${\dot{M}}_{g}$ does plays a relatively minor role in Equation (9), resulting in a tight relationship between SFR and ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$. Although it is important to bear in mind that we have made a simple estimate of ${\epsilon }_{\mathrm{MD}}$ which may not be the actual mass deposition efficiency, our analysis demonstrates how the relationship between ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ and star formation may be used to constrain the processes governing the cooling and condensation of gas in the ICM.

In Fogarty et al. (2015), we attempted to estimate the cooling rate of the ICM by measuring ${\dot{M}}_{g}$ for gas that was at a radius $\lt 35\,\mathrm{kpc}$ or for gas that had an average ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ ratio below 70.¹⁰ We found that while reddening-corrected UV photometric SFRs scale with ${\dot{M}}_{g}$ in both cases, star formation appears to be increasingly "inefficient" as the SFR decreases. The lowest SFRs we observed were $\sim 0.1 \% \mbox{--}1 \%$ of ${\dot{M}}_{g}$, while at the other extreme the two quantities were comparable. Such a trend would be expected if, as Equation (9) indicates, ${\epsilon }_{\mathrm{MD}}$ scales with the SFR.

5.1.1. The Role of t_ff in the SFR–${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ Relationship

Recent work examining the critical condition for the onset of BCG activity in a cool-core cluster suggests that this activity is driven by t_cool, not ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ (Hogan et al. 2017). These results may imply that we ought to observe a relationship between SFR–t_cool without any significant contribution from t_ff. An SFR–t_cool relationship may be observed if the dominant driver of condensation is the AGN jet propelling low-cooling-time gas from lower to higher altitudes, where it can condense. Since the cooling time of the uplifted gas determines the radius where the gas becomes thermally unstable, it also determines how much work must be done by the jet to lift the gas to a radius where condensation can occur. Therefore, if the uplift of plasma is driving condensation, we expect to see an SFR–mass condensation efficiency relationship in the form of an SFR–t_cool relationship. Such a relationship would appear similar to an SFR–${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ relationship in a sample such as CLASH with a narrow range of cluster masses, since all of the clusters in CLASH have similar freefall times in their cores.

We thus investigated the possibility that the underlying relationship we observe in Figure 5 is primarily between SFR and t_cool. The range of t_ff at $0.025\,{R}_{500}$ in CLASH clusters with star-forming BCGs is smaller than the range of t_cool—while t_ff varies by a factor of ∼1.5 across the sample, t_cool varies by a factor of ∼3. Hence, the effect of t_ff on the relationship with SFR is modest, and our ability to distinguish between an underlying SFR–${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ relationship versus an SFR–t_cool relationship is limited. Several lines of reasoning provide evidence for a non-negligible contribution from t_ff, although we cannot rule out the interpretation that the underlying relationship we observe is solely between SFR and t_cool.

SFR and t_cool are correlated as shown in Figure 9, with a Spearman correlation coefficient of −0.90 (and a Pearson coefficient of −0.90). We calculate an intrinsic scatter between the two quantities of ${0.05}_{-0.01}^{+0.03}$ dex ($\lt 0.18$ dex at $3\sigma$). The Spearman correlation coefficient in the SFR–${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ and SFR–t_cool relationships as a function of sampling radius are presented in Figure 10. We do not find a substantial difference between t_cool and ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ in terms of how tightly these quantities relate to the SFR. However, the distinction between t_cool in the star-forming and non-star-forming CLASH clusters is less clear than the distinction between ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ in these two populations, and t_cool for two clusters with non-star-forming BCGs is comparable to t_cool in the clusters with BCGs exhibiting $0.1\mbox{--}1\,{M}_{\odot }$ yr⁻¹ of star formation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The situation becomes clearer when we examine the relationship between ${t}_{\mathrm{cool}},{t}_{\mathrm{ff}}$, and the residuals in the fits to the CLASH SFR–${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ and SFR–t_cool data sets. In Figure 11, we show t_ff versus the residuals when we fit a log–log relationship to the SFR–t_cool data set at a radius of $0.025\,{R}_{500}$ and at $0.075\,{R}_{500}$. We also show t_ff versus the residuals of SFR–${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ at $0.025\,{R}_{500}$ and t_cool versus the residuals of SFR–${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ at $0.025\,{R}_{500}$. We find that at $0.025\,{R}_{500},{t}_{\mathrm{ff}}/\langle {t}_{\mathrm{ff}}\rangle$, where $\langle {t}_{\mathrm{ff}}\rangle$ is the mean t_ff for the sample, is correlated with the residuals in the fit to the SFR–t_cool data set, $\tfrac{{t}_{\mathrm{cool}}}{\mathrm{Predicted}\ {t}_{\mathrm{cool}}}$, where ${\text{}}{Predicted}\ {t}_{\mathrm{cool}}$ is the t_cool predicted by the SFR–t_cool relationship for a given SFR. These two quantities are consistent with $\tfrac{{t}_{\mathrm{cool}}}{\mathrm{Predicted}\ {t}_{\mathrm{cool}}}={t}_{\mathrm{ff}}/\langle {t}_{\mathrm{ff}}\rangle$, implying that dividing t_cool by t_ff will offset the residuals in the SFR–t_cool relationship. The relationship between t_ff and $\tfrac{{t}_{\mathrm{cool}}}{\mathrm{Predicted}\ {t}_{\mathrm{cool}}}$ at $0.025\,{R}_{500}$ has a positive slope with $\sim 98 \%$ confidence. The plot in Figure 11 does not have obvious outliers, so the scatter reduction seen at $0.025\,{R}_{500}$ by dividing t_cool by t_ff is not attributable to reducing the residuals in an extreme outlier. There is also no evidence that t_ff or t_cool are correlated with the other residuals we examined.

Zoom In Zoom Out Reset image size

We find a relationship between ${\rm{\Delta }}{t}_{b}$ and t_cool that shows a possible positive correlation between these two quantities with ${\rm{\Delta }}{t}_{b}$ approaching t_cool at gigayear timescales. Since ${\rm{\Delta }}{t}_{b}$ is potentially correlated with the SFR, and individual measurement uncertainties are relatively large (∼0.3–0.4 dex) compared to the range of ${\rm{\Delta }}{t}_{b}$ values measured in the sample (1.5 dex), the observed relationship has limited power to constrain models of cluster-scale feedback evolution. However, our measurements of the starburst durations have implications for understanding how AGN-regulated feedback progresses over time.

Four of the eleven BCG starburst durations are reported as lower limits. These BCGs may be undergoing continuous star formation. However, given the association between longer ${\rm{\Delta }}{t}_{b}$ and lower SFR, and shorter ${\rm{\Delta }}{t}_{b}$ and higher SFR, we suspect that star formation in these BCGs is slowly decaying. Alternatively, BCGs undergoing feedback may initially exhibit large starbursts before settling down to a relatively steady state with anywhere between ∼1 and 10 ${M}_{\odot }$ yr⁻¹ of star formation. Star formation may also "flicker" on timescales that are short relative to the values of ${\rm{\Delta }}{t}_{b}$ we measure.

Along with recent spectroscopic observations in the far-UV, the long burst durations we observe suggest that star formation and the thermodynamical state of the ICM in the cool core is temporally decoupled from the gas condensation rate fueling AGN feedback, which may feature irregular spikes and dips over time. Recent observations of RXJ 1532.9+3021 with the Cosmic Origins Spectrograph (COS; Green et al. 2012) reveal that star formation in this BCG exceeds the upper limit on gas cooling measured by N v and O vi by a factor of ∼10, suggesting that the initial buildup of molecular gas in this system has run its course (Donahue et al. 2017). Meanwhile, similar spectroscopy in Abell 1795 and the Phoenix cluster reveal that gas cooling outstrips the SFR (McDonald et al. 2014a, 2015). Delayed consumption of molecular gas may allow the star formation history of BCGs to "smooth over" these intermittent spikes. Slowly decaying BCG starbursts may also trace a transition from an initial mode of rapid condensation (with ∼100–1000 ${M}_{\odot }$ yr⁻¹ of gas condensing) to a long-duration mode of more modest condensation (with ∼1–10 ${M}_{\odot }$ yr⁻¹ of condensation). If star formation lags behind cooling in making this transition, the resulting observables would be consistent both with our results and the offset between SFRs and cooling rates seen in Molendi et al. (2016) and Donahue et al. (2017).

A positive correlation between ${\rm{\Delta }}{t}_{b}$ and t_cool and the existence of long-lived ($\gt 1$ Gyr) starbursts are consistent with the mass condensation efficiency interpretation of the ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$–SFR scaling relationship discussed in Section 5.1. Assuming that the duration of the BCG starburst is a proxy for the duration of AGN feedback, then longer periods of star formation in a BCG imply more energy injection into the surrounding ICM, which in turn raises the cooling time of the surrounding ICM. This in turn diminishes the rate at which molecular gas condenses, and gradually shuts off both the starburst and feedback. The duration ${\rm{\Delta }}{t}_{b}$ and t_cool ought to converge in the limit of a long starburst, although feedback mechanisms could run out of fuel and shut down before the two quantities converge depending on how inefficiently mass condenses out of the ICM. Alternatively, a quasi-steady state may be reached with effectively continuous star formation in the BCG coupled with ICM plasma with t_cool of a few gigayears being replenished by hotter gas at about the same rate it condenses. In this case, star formation and feedback may decay very slowly, or not at all. Our results do not rule this latter scenario out, leaving the following possible options: (1) BCG activity is cyclical, (2) BCG activity approaches a slowly decaying quasi-steady state, or (3) BCG activity is cyclical but the duty cycle is comparable to the age of the cool core of the cluster.

Finally, it has been noted that there are few examples of BCGs with post-starburst spectra in most BCG samples (Liu et al. 2012; Loubser et al. 2016). Indeed, none of the SOAR or SDSS spectra of CLASH BCGs in Fogarty et al. (2015) show post-starburst features. The lack of post-starburst BCGs is consistent with our estimates of ∼Gyr duration episodes of star formation. Wild et al. (2009) notes that for a galaxy with an exponentially decaying starburst to feature a post-starburst spectrum, the burst's decay timescale would have to be $\lesssim 0.1\,\mathrm{Gyr}$ and the burst would have to account for at least 5%–10% of the galaxy's stellar mass. Given the nature of long-duration star formation in our observations and models of feedback-regulated cooling, we would not expect to see a substantial population of post-starburst BCGs.

Large BCG starbursts (with SFRs $\gtrsim \,100\,{M}_{\odot }$ yr⁻¹) may differ from their more modest counterparts in several ways. First, the SFRs and dust masses in the CLASH BCGs are consistent with using Hjorth et al. (2014) to describe the relationship between BCG star formation and dust mass. Specifically, MACS 1931.8−2653 and the Phoenix cluster match the flattening slope and eventual turnover in the Hjorth et al. (2014) relation closely, raising the possibility that these BCGs harbor starbursts similar to the starbursts in the submillimeter-detected population of galaxies studied by Hjorth et al. (2014). The Hjorth et al. sample extends the study of dust and star formation conducted in da Cunha et al. (2010) to cover starbursts to galaxies forming stars at $\gtrsim 1000\,{M}_{\odot }$ yr⁻¹. Since the largest BCG starbursts in the CLASH sample also began forming stars more recently, we hypothesize that, like these massive field galaxy starbursts, BCGs with large star formation rates are either forming or building up their dust reservoirs.

Second, both MACS 1931.8−2653 and RXJ 1532.9+3021 boast prodigious SFRs and distinctive X-ray cavities. Furthermore, MACS 1931.8−2653 exhibits an X-ray-loud AGN, a feature that is not obvious in the other CLASH clusters. These features are noteworthy because the starbursts in RXJ 1532.9+3021 and MACS 1931.8−2653 are the youngest in our sample at ${\mathrm{log}}_{10}{\rm{\Delta }}{t}_{b}=-{0.39}_{-0.41}^{+0.38}$ and $-{1.01}_{-0.34}^{+0.35}$ Gyr, respectively. Our findings suggest that stronger X-ray features may be associated with younger, larger BCG starbursts, such as the extreme example of BCG star formation present in the Phoenix Cluster (McDonald et al. 2012, 2013, 2014b).