The Evolution of AGN Activity in Brightest Cluster Galaxies

Galaxy clusters are the most massive gravitationally bound and collapsed objects in the universe (Voit 2005). Because of their extremely deep potential wells, the temperature of the intracluster medium (ICM) is high enough to radiate X-rays. The central parts of clusters, which have the densest X-ray emitting gas, often have shorter cooling times than the Hubble time, implying that the hot X-ray gas should have had enough time to cool and form large inward flows of cooling material, known as cooling flows (Sarazin 1986; Fabian 1994). However, multiwavelength observations have only seen a fraction of the massive cooling flows that are expected from standard cooling models (e.g., O'Dea et al. 2008; Donahue et al. 2015; McDonald et al. 2018). This is referred to as "the cooling-flow problem," and active galactic nucleus (AGN) feedback is thought to be responsible for preventing the hot gas from cooling by propagating energy from the supermassive black hole (SMBH) to the ICM. The two primary modes of AGN feedback are the kinetic mode, in which relativistic jets push the hot gas aside and create cavities, and the quasar mode, in which the radiation comes from the accretion disk (see reviews by Fabian 2012; McNamara & Nulsen 2012).

With the recent development of galaxy cluster surveys that use the Sunyaev–Zel'dovich (SZ) effect (Sunyaev & Zeldovich 1972), such as the South Pole Telescope (SPT; Carlstrom et al. 2011; Bleem et al. 2015, 2020; Huang et al. 2020) and the Atacama Cosmology Telescope (ACT; Hilton et al. 2018, 2021), the number of known high-z (z > 1) clusters with good mass estimates has increased dramatically. This has enabled many studies of the evolution of AGN feedback in clusters over cosmic time (McDonald et al. 2013, 2017; Gupta et al. 2020). However, the evolution of AGN feedback in galaxy clusters with redshift remains poorly understood. In particular, Hlavacek-Larrondo et al. (2015) found no evidence for evolution in jetted power generated by AGN feedback from X-ray cavities over the past 7 Gyr (z = 0.8). An earlier study by Hlavacek-Larrondo et al. (2013) suggested that the fraction of brightest cluster galaxies (BCGs) with X-ray bright nuclei is decreasing with time (or increasing with redshift), suggesting a strong evolution in radiative mode feedback. In contrast, a recent study looking for nuclear BCG X-ray emission in Chandra archival data instead found no evidence for an evolution between two redshift bins (〈z〉 ∼ 0.25 and 〈z〉 ∼ 0.65; Yang et al. 2018). The disagreement between various studies about the evolution of AGN feedback restricts our ability to fully understand this issue.

In this work, we calculate the AGN-hosting BCG fraction by identifying BCGs in the 2500 deg² SPT-SZ cluster samples (Bleem et al. 2015) and by classifying whether they are AGNs based on mid-IR data. The SZ cluster catalogs allow for an effectively mass-selected sample of clusters, making it possible to study the evolution of galaxy clusters over time. In addition, the SZ catalogs typically have less contamination compared to optical/IR catalogs. The fraction of BCGs hosting luminous AGNs is an important indicator for AGN fueling processes, availability of cold clumpy gas in the centers of clusters, the duration and duty cycle of the AGNs, and how BCGs and the host clusters grow and evolve together. This is because additional physical mechanisms are often required to explain the transport of the cold gas, which serves as the primary fuel source for the central black holes. The fact that we find a relative absence of AGNs in the centers of clusters has led us to study many physical processes, including ram pressure stripping (Gunn & Gott 1972), tidal effects from the cluster gravitational potential (Merritt 1983), and the lack of new infall of cold gas (Larson et al. 1980). Similarly, understanding the evolution of AGN activities in BCGs will help us understand the evolution mechanism of galaxy clusters, and how the feedback might play a role in it.

Our goal for this paper is to study the redshift evolution of the AGN-hosting BCG fraction up to z = 1.3 to understand the fueling processes in the centers of clusters, determine when AGN feedback is fully established, and identify whether there are any more extreme AGN-hosting BCGs in the sample, similar to the Phoenix cluster. The paper is organized as follows. In Sections 2 and 3 we summarize the data and additional information used in this paper. The results and their implications are presented in Sections 4 and 5, respectively. We conclude our work in Section 6. We assume H₀ =70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_λ = 0.7. All errors are 1σ unless noted otherwise.

2.1. The SPT-SZ 2500 deg² Cluster Sample

2.2. Position of BCGs

Given the diversity of BCG colors, morphologies, and assembly state as a function of redshift, typical identification algorithms may be biased when they select BCGs based on single-band fluxes. We have instead developed a novel BCG identification pipeline that uses the full probability distribution of redshift and stellar mass for every object within 500 projected kiloparsec of the SZ cluster center to assign BCG likelihoods. Photometry is provided by the Dark Energy Survey (Year 3) catalogs (Jarvis et al. 2021; Sevilla-Noarbe et al. 2021), cross-correlated with unWISE (Lang 2014), which is a combination of WISE and NEOWISE images. Various cuts and flags are used to avoid stars and objects with poor photometric measurements. WISE is an IR-wavelength space satellite with four IR filters, including W1 (λ_cen = 3.6 μm), W2 (4.3 μm), W3 (12 μm), and W4 (22 μm; Wright et al. 2010). The satellite operated for two years with cryogen until 2011, before it was reactivated and resumed operations as NEOWISE in 2013. It has continued to observe ever since (Mainzer et al. 2014).

Probability distributions of photometric redshift and stellar mass for each source are estimated with EAZY (Brammer et al. 2008) and FAST (Kriek et al. 2009), respectively. We then randomly sample from each distribution to find the most massive cluster galaxy at the cluster redshift within each field, iterating this process 10⁵ times to build up a BCG likelihood for each galaxy. In this way, all galaxies are assigned a value between 0% and 100% probability of being the BCG within each cluster. Full details about the pipeline, along with the BCG catalog, will be provided in Noble et al. (in preparation). The top three panels of Figure 1 show optical images of example SPT galaxy clusters with identified BCGs in white squares. This demonstrates that the algorithm selects likely BCGs that match the galaxies that typical or traditional visual BCG identification methods would select over a wide range of redshift.

Zoom In Zoom Out Reset image size

2.3. Data for AGN Selection

Most photometric techniques for identifying AGNs are severely biased toward unobscured (type 1) AGNs since their nuclear emissions dominate the emission of the host galaxies, making these AGNs easily identifiable. This implies that most obscured (type 2) AGNs are often underrepresented in most studies. The most promising techniques for identifying both obscured and unobscured AGNs include radio, hard X-ray, and mid-IR selections. However, not all AGNs are radio loud (e.g., Stern et al. 2000), and the current hard X-ray satellites remain limited in their sensitivity and field of view. This leaves mid-IR selection as a popular technique to quickly identify large AGN populations (obscured and unobscured). The idea of mid-IR selection is to separate between the power-law AGN spectrum and the blackbody stellar spectrum of galaxies, which peaks at rest frame 1.6 μm. The power-law spectra of the AGNs are due to the thermal emission from the warm-hot dust in the torus, which is heated by absorbing shorter-wavelength photons from the accretion disk (Stern et al. 2012; Hickox & Alexander 2018). This implies that the emission is not strongly suppressed by the dusty torus, unlike at the UV to near-IR wavelength, allowing this technique to detect more obscured AGNs. Additionally, with the first all-sky data release of the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) in 2012, mid-IR selection became one of the top methods of probing the AGN population over the entire sky without additional observations.

One drawback of the mid-IR selection technique is that the host galaxy is still bright at these wavelengths, limiting the detection of low-luminosity AGN, which have to compete with a bright stellar continuum (Stern et al. 2012). This means that AGNs that are selected by their mid-IR color tend to be brighter relative to the host galaxies than those selected by other techniques. For example, Assef et al. (2013) found that in the sample of relatively luminous AGNs (L_AGN/L_host > 0.5), the luminosity of mid-IR AGNs tends to be greater than ∼5 × 10⁴⁴ erg s⁻¹, taking into account the bolometric correction from Singal et al. (2016). Assuming the efficiency of turning accreting matter into energy _acc = 0.1 and a typical mass of the SMBH M_BH ∼ 10⁸–10⁹ M_⊙ (Russell et al. 2013), the black hole mean accretion rate (${\dot{M}}_{{BH}}={L}_{\mathrm{bol},\mathrm{nuc}}/{\epsilon }_{\mathrm{acc}}{c}^{2}$) of mid-IR selected AGNs should be greater than 4 × 10⁻³–$4\times {10}^{-2}\,{\dot{M}}_{\mathrm{Edd}}$, where ${\dot{M}}_{\mathrm{Edd}}$ is the limiting Eddington accretion rate. This level of accretion is relatively high compared to typical optical/radio AGNs, which have an accretion rate of about 10⁻⁶–${10}^{-2}\,{\dot{M}}_{\mathrm{Edd}}$ (McDonald et al. 2021). This implies that the mid-IR technique will mainly identify a brighter and more massive AGN population. From now on, AGNs mentioned in this paper mean the mid-IR selected AGN population.

2.3.1. Mid-infrared Data from WISE

Instead of using the main source catalog from WISE (AllWISE), which only includes the data obtained from 2010 to 2011, we make use of CatWISE to obtain the best photometry with available data. CatWISE is an updated all-sky IR source catalog that combines the 2010 and 2011 data from WISE with the 2013 through 2016 NEOWISE data (Eisenhardt et al. 2020). The caveat is that the CatWISE catalog only includes 3.6 (W1) and 4.3 (W2) μm data. In this work, we use the Preliminary CatWISE catalog to obtain the WISE color for each BCG, which has the advantage that more than four times as many exposures are included as were used for the AllWISE catalog while using the same AllWISE software. This makes a comparison with previous works more straightforward.

For every identified BCG (with a probability >5%) from Section 2.2, we search for mid-IR counterparts in the CatWISE catalog within a radius of 3'' from the identified BCG because the typical FWHMs for W1 and W2 are 608 and 684, respectively. Both W1 and W2 are converted from Vega into AB magnitudes using the correction from the Explanatory Supplement to the WISE Products. ⁵⁷ The bottom three panels of Figure 1 show W1 images of the three SPT clusters, showing that their BCG candidates can be detected with WISE data.

2.3.2. CatWISE Color Correction

We perform a comparison test between the AllWISE and CatWISE catalogs. The test is carried out by comparing the (W1 − W2)_AB color of bright objects (16.5 < W1_AB < 18) between the two catalogs within 10 arcmin of one field. A histogram of the color differences is plotted in Figure 2, showing the offset between the color from AllWISE and CatWISE to be around 0.042. Further investigation reveals that this is due to the gradual diminishing of the W2 throughput with time, leading to a bluer (W1 − W2)_AB color compared to AllWISE. We apply this correction of 0.042 mag to CatWISE W1–W2 colors.

Zoom In Zoom Out Reset image size

3.1. AGN Selection

With the WISE satellite, Stern et al. (2012) developed a well-known formula for quickly identifying AGN candidates with a simple color criterion, (W1 − W2)_Vega ≥ 0.8 or (W1 − W2)_AB ≥ 0.16. One benefit of mid-IR selected samples is that both unobscured (type 1) and obscured (type 2) AGN can be identified (Lacy et al. 2004; Stern et al. 2005, 2012). However, because the colors of galaxies drastically change over a large redshift range, this simple criterion is not accurate enough to characterize a large population of AGNs. To increase the number of AGNs we can identify, we develop a new color criterion that depends on the redshifts of galaxies.

3.1.1. EzGal Galaxy Color Model

To determine whether each BCG harbors an AGN, we calculate the expected color for typical elliptical galaxies using EzGal. ⁵⁸ EzGal calculates the magnitude evolution as a function of redshift from evolving the spectral energy distribution (SED) models of a stellar population with time and projecting them through filters (Mancone & Gonzalez 2012). This calculation takes into account both the stellar evolution of a galaxy as young stars evolve and the wavelength shift due to the distance of a galaxy. To find the model that best describes our overall sample, we perform a grid search between three stellar population model sets (i.e., Bruzual & Charlot 2003; Maraston 2005; Conroy et al. 2009), various formation redshifts, two different initial mass functions (IMF) (i.e., Salpeter 1955; Chabrier 2003), star formation history as a single exponential decaying burst of star formation with an e-folding time parameter (τ) between 0.1 and 10 Gyr, and the representative metallicity (Z) for our galaxy sample from 0.001 to 0.03. Ultimately, the best-fit model (based on the chi-square test) is a Bruzual & Charlot (2003) stellar model with a formation redshift of (z_f ) 3.5, a Salpeter (1955) IMF, τ = 0.1 Gyr for the star formation history, and a metallicity of 0.016. The solid orange line in the top panel of Figure 3 shows the expected W1–W2 color evolution, generated from the EzGal model with this particular set of parameters. The bottom panel shows the residual from the expected value of W1-W2 for each BCG. It demonstrates that the scatter is distributed around zero with a relatively weak redshift dependence, implying that we have successfully removed the continuum contribution. In this work, galaxies that are redder (a residual greater than 0.2) than typical elliptical galaxies based on the EzGal model are considered AGN candidates. For a range of threshold values from 0.2 to 0.4, the highest redshift bin has more AGN-hosting BCG sthan the lowest redshift bin, implying that an increase in the fraction of AGN-hosting BCGs is independent of this choice. To further test this notion, we perform a heteroscedasticity test, specifically, the Breusch-Pagan test, in the bottom panel of Figure 3, which shows whether the scatter of the IR residual depends on redshift, regardless of the choice of a threshold value. The test results in a P-value = 0.0047, meaning that the heteroscedasticity is present and the scatter of W1–W2 residuals depends on the redshift, implying that an increase in the fraction is feasible regardless of the threshold choice.

Zoom In Zoom Out Reset image size

One assumption that we apply in this section is that we only consider a single-burst stellar population model with a single formation redshift, star formation history, and metallicity. In Figure 4 we consider both a single-burst stellar population model and a more complicated two-age stellar population (old and young) model with a wide range of parameters for both models. The model with two stellar populations keeps the same parameter sets from a single-burst model for the "old" population, while a "young" population is represented by a 50 Myr old stellar population at all redshifts. Even though these two models are likely not sufficient to describe our data, more sophisticated models would be unconstrained by the data that we have available. Based on these single-age and two-age models, we find no combination of formation time, metallicity, and IMF that can fully account for the observed evolution in the mid-IR excess, as shown in the right panel of Figure 4. We propose that this mid-IR excess comes from a dusty torus, which is a signature of an actively accreting SMBH. It is difficult to imagine other astronomical sources for this emission because star formation typically yields significantly cooler dust temperatures with a peak brightness of ∼100 μm instead of ∼1–10 μm. On the other hand, it could also be that our current population models are not adequate to describe the data. We should keep this caveat in mind when we discuss the implication of our results.

Zoom In Zoom Out Reset image size

As a test to see how a starburst can affect the mid-IR color, we consider Abell 1835 (Ehlert et al. 2011) and RX J1532.9+3021 (Hlavacek-Larrondo et al. 2013), which are the most star-forming BCGs known (SFR ∼ 100 M_⊙ yr⁻¹; McDonald et al. 2018) that also lack evidence of a strong AGN. The two colored stars in Figure 3 demonstrate that even though a star-forming BCG would have boosted mid-IR emission due to dust, polycyclic aromatic hydrocarbon molecules (PAH), and molecular gas, the emission is not as strong as the power-law spectra of AGNs, and our selection does not include these two BCGs. On the other hand, the two clusters with the most luminous AGNs (H 1821+643; Russell et al. 2010, and IRAS 09104+4109; O'Sullivan et al. 2012) are easily detected with our criterion.

3.1.2. Spitzer Color Verification

Because the point-spread functions (PSFs) of the two WISE bands are not small (PSF _W1 = 608 and PSF _W2 = 684), we compare the results from WISE mid-IR color with those from the Spitzer Space Telescope. Spitzer is an IR telescope with the Infrared Array Camera (IRAC; Fazio et al. 2004) as one of its main science instruments. IRAC is a four-channel imaging camera capable of taking simultaneous images at wavelengths of 3.6, 4.5, 5.8, and 8.0 μm. Thus, channels 1 and 2 ([3.6] and [4.5]) on IRAC are roughly equivalent with W1 and W2 from WISE, but with the benefit of having a much better PSF at 195 and 202, respectively (Fazio et al. 2004).

A certain fraction (∼35%, predominantly at z > 0.8) of the SPT cluster sample has been observed with IRAC. For verification, we compare the [3.6]–[4.5] colors of our AGN-hosting BCG candidates with their W1–W2 colors. If the Spitzer color, which has a higher angular resolution, is bluer (smaller) than the WISE color, it shows that there is a contamination from nearby galaxies within the WISE aperture. On the other hand, if the Spitzer color is redder (larger), it implies that the object is even more likely to be an AGN. Figure 5 shows the comparison between WISE's W1–W2 (gray squares) and Spitzer's [3.6]–[4.5] color (circles) for our AGN candidates. We find that most AGN candidates have a difference of Spitzer and WISE color that is either compatible (60% of the sources have a difference of within ±0.07 mag, which is roughly the mean of the WISE color uncertainty) or that Spitzer is slightly redder (33% of the sources have a Spitzer color larger than the WISE color by ∼0.2 mag). This suggests that most of our AGN candidates are likely to be real quasars. One clear exception is SPT-CL J2146-4633, which has a WISE color much redder than Spitzer. Further investigation shows that there is a point-like source near the location of the object, but not at the BCG location, meaning that the WISE color is probably contaminated by a nearby AGN, while the Spitzer color is not. This object has been removed from the further analysis.

Zoom In Zoom Out Reset image size

3.2. Calculating the AGN fraction

Table 1 shows the number of AGN-hosting BCG in the four redshift bins, whereas the blue points in Figure 6 show the fraction of AGN-hosting BCGs in the four bins with their corresponding 68% confidence intervals. We observe that the fraction is increasing with redshift in the SPT sample. Gray crosses show the fraction of points that have residuals lower than -0.2. If the scatter is truly random, there should be a similar number of points below -0.2 compared to those above +0.2, which are classified as AGNs. This result implies that the increasing trend of the AGN-hosting BCGs fraction is not a result of the data quality. Such a trend has been suggested and shown in previous works (e.g., Hlavacek-Larrondo et al. 2013; McDonald et al. 2016; Bîrzan et al. 2017; Gupta et al. 2020; Mo et al. 2020). In particular, we show that the fraction is ∼2% at z ≈ 0.5, which is consistent with what Somboonpanyakul et al. (2021a) found from studying extreme central BCGs in clusters. We note that because some AGN-dominated galaxies will have poor photometric redshift constraints as they are estimated from the stellar spectrum and not from the power-law spectrum of the AGN, we might misidentify these galaxies in our BCG-finding algorithm. This implies that the number of BCGs with central AGNs found in this work gives a lower limit for the AGN-hosting BCG fraction, and the actual evolution could be even stronger.

Zoom In Zoom Out Reset image size

One way to evaluate the observed result with the chosen bins is to perform Fisher's Exact Test in order to see whether there are any associations between two categorical variables. The result shows that we can reject the null hypothesis of independence with P-value = 0.00045, meaning that there is a statistically significant association (>99.9%) between redshift and whether or not BCGs host AGNs.

The approaches taken in this work are as follows: (i) we include all BCG candidates with a probability higher than 20% in our sample instead of picking only one BCG per cluster, and (ii) we calculate the fraction of BCGs with AGNs over the total number of BCG candidates, and not the fraction of clusters with AGN-hosting BCGs over the total number of clusters. The reason for these two assumptions is that we want to include AGN-hosting BCGs from systems with more than one obvious BCG, which are typical for merging systems such as the Coma cluster (Zwicky 1933), and the Bullet cluster (Markevitch et al. 2004). We perform consistency checks to address both of these assumptions. Figure 7 shows the fraction of AGN-hosting BCGs when we consider the most likely BCG candidates, every BCG candidate with a probability higher than 20%, and every BCG candidate with a probability higher than 10%. This figure shows that the increasing trend of the fraction of AGN-hosting BCGs over redshift remains consistent in all three scenarios, regardless of how we select BCGs.

Zoom In Zoom Out Reset image size

On the other hand, Figure 8 shows the results when we use different definitions of AGN fractions. The gray and blue points in Figure 8 are calculated with the total number of clusters as a denominator instead of the number of BCG candidates. For the blue points, we consider one BCG per cluster and include both the probability of being BCGs, as calculated in Section 2.2, and the uncertainty of the mid-IR color for each BCG to emphasize the fact that the uncertainties of identifying BCGs and BCG colors are higher at high redshift. The empty gray dots are the largest possible fractions, which are calculated from clusters for which any of their potential BCGs are considered as AGNs, while the empty gray squares are the smallest possible fractions by counting only clusters for which all of their BCG candidates are classified as AGNs. This figure illustrates that all of these definitions qualitatively give the same conclusion as our initial results.

Zoom In Zoom Out Reset image size

We compare our results with the AGN fraction in field galaxies to determine whether there is a difference in the fractions between the two environments. The green and pink squares in the left panel of Figure 9 show the field X-ray AGN fractions from the zCOSMOS survey (Silverman et al. 2009) and the Chandra Multiwavelength Project results (ChaMP; Haggard et al. 2010), respectively. The results from our work are consistent with these two results, suggesting that the source of fuel required for AGN accretion in field galaxies could be similar to that in the brightest cluster galaxies. Additional evidence for the AGN fraction evolution in field galaxies has been seen in other works. For example, Lehmer et al. (2007) reported an evolution in early-type galaxies (z ∼ 0.7) consistent with the (1 + z)³ pure luminosity evolution model. The dotted gray line in Figure 9 shows the curve for (1 + z)³, although it is only intended to be illustrative because it is arbitrarily normalized. The dashed line instead shows the curve for (1 + z)^5.3. This line was first suggested by Martini et al. (2009), who showed that the AGN fraction of cluster members increases as ∼(1 + z)^5.3 for AGN with an X-ray luminosity higher than L_x > 10⁴³ erg s⁻¹, hosted by luminous galaxies. We also fit the power-law model (∝(1 + z)^α ) to the blue points in Figure 9 and find a power-law exponent α = 4.1 ± 1.0, as shown by the dash-dotted brwon line, which is consistent with the results from both Lehmer et al. (2007) and Martini et al. (2009). Nevertheless, there are caveats regarding the relation between cluster BCGs and field galaxies. One concern is that the AGN selection criteria for both BCGs and field galaxies are different, making it difficult to make a direct comparison between the two. In addition, according to the work about the evolution of AGN luminosity, which shows that AGNs in galaxies tend to be brighter at high redshift (Hasinger et al. 2005; Silverman et al. 2008), we would naturally expect to find a higher AGN fraction at high redshift because we usually selected AGNs based on a certain luminosity threshold.

Zoom In Zoom Out Reset image size

The results obtained in Section 4 demonstrate that given a single-burst stellar population model, there is an increase in the fraction of AGN-hosting BCGs with redshift. This finding is consistent with previously published studies (e.g., Hlavacek-Larrondo et al. 2013; Bîrzan et al. 2017), which focused on different samples with distinct selection effects. In particular, the results from Hlavacek-Larrondo et al. (2013) and ours, as shown in the middle panel of Figure 9, show the same trend of an increasing fraction of AGN-hosting BCGs with redshift. However, the normalizations are vastly different. Hlavacek-Larrondo et al. (2013) claimed that the fraction of active BCGs is 30% at z ≈ 0.1 and 60% at z ≈ 0.5. On the other hand, we show that the fraction is smaller than 20% at all redshift bins. A possible explanation stems from the very different samples that these works consider. While this study is based on an effectively mass-selected sample of clusters, without any consideration of X-ray properties, the sample used by Hlavacek-Larrondo et al. (2013) focused solely on highly X-ray luminous clusters (L_X,cluster > 3 × 10⁴⁴ erg s⁻¹), which show clear X-ray cavities. To place these two studies on the same scale, we modify the denominator used by Hlavacek-Larrondo et al. (2013) in calculating the AGN fraction to account for the full parent population of clusters from which their sample of 32 clusters was drawn. Given that this previous work included a subsample of clusters drawn from the REFLEX (Böhringer et al. 2004), eBCS (Ebeling et al. 2000), MACS (Ebeling et al. 2001), and SPT-XVP (McDonald et al. 2013) surveys, we consider these surveys in their entirety as the total population (the denominator) when calculating the AGN fraction. With this rescaling, we find that the fraction of BCGs hosting powerful AGN is consistent between our work and Hlavacek-Larrondo et al. (2013), and that the observed evolution is consistent in the two studies. We note that AGN are identified in very different ways in these two samples, but they both select at the extreme end—the few most mechanically powerful and the most IR-bright outbursts. The fact that these two works agree after the aforementioned rescaling is applied is reassuring, although a much more thorough analysis (and proper bias correction) is needed before conclusions about the coevolution of jet power and mid-IR emission can be made.

An increase in the fraction of BCGs hosting central AGNs with redshift suggests that the accretion rates of the SMBHs in the BCGs are higher at high redshift because the AGN luminosity is proportional to the accretion rate. Several works about the relation between the mean black hole accretion rates and the cavity (kinetic)/quasar (radiative) power of the central AGN (Churazov et al. 2005; Russell et al. 2013) have shown that as black hole accretion increases in the BCGs, the cavity power of the AGN also increases to counteract the cooling from the accretion in the form of a negative feedback cycle. However, as the black hole accretion rates rise to near the Eddington limit, the cavity/jet power seems to be saturated, and the radiative power tends to dominate at this level of accretion. The fact that the radiative power from AGNs usually promotes more cooling in the ICM instead of preventing it suggests that a well-regulated feedback system between the central black hole and its host cluster is no longer possible at high accretion, implying that some galaxies might not have a fully established AGN feedback loop at this redshift range.

A similar conclusion has been reached in the work related to star-forming galaxies (Webb et al. 2015; McDonald et al. 2016; Bonaventura et al. 2017), which showed that the fraction of starburst BCGs is higher at high redshift (z > 1). Specifically, McDonald et al. (2016) found the fraction of BCGs with SFR over 10 M_⊙ yr⁻¹ to be 34 ± 5% at 0.25 < z < 1.25, compared to ∼1%–5% at z ∼ 0. The right panel of Figure 9 compares the fraction of AGN-hosting BCGs in this work with the fraction of starburst BCGs (SFR_BCG > 100 M_⊙ yr⁻¹) from McDonald et al. (2016), demonstrating that in the center of cluster environments, both massive starburst galaxies and bright AGNs behave similarly. These two results strongly indicate that AGN feedback might not be as effective to prevent overcooling at high redshift as we have previously thought.

All of these results lead us to suspect that the reason for the observed redshift trend and the breakdown of AGN feedback at high redshift comes from the fact that there is an abundance of cold gas at that redshift. In the local universe, BCGs typically grow by merging with gas-poor satellites without triggering any AGN activity. However, BCGs at high redshift could grow by merging with gas-rich members instead. Cold gas from the mergers could be a source of fuel for increasing AGN activities in the center of clusters. This is consistent with the picture we obtain from the studies of starburst BCGs (Webb et al. 2015; McDonald et al. 2016) because cold dense clouds from gas-rich mergers could provide enough matter required for creating new stars. Further evidence supporting a gas-rich merger explanation includes the prevalence of cluster galaxies with massive CO or cold gas reservoirs at high redshift (Noble et al. 2017, 2019; Hayashi et al. 2018; Markov et al. 2020) and the detections of molecular gas in many BCGs (Dunne et al. 2021). This scenario can also explain recent studies about the cool-core (CC) fraction, which shows no sign of evolution over the same redshift range (McDonald et al. 2017; Ruppin et al. 2021). If AGN feedback breaks down at high redshift, one would expect that the CC fractions of clusters would be higher because more gas should have been cooled near the center. However, if black hole accretion and star formation in high-redshift BCGs are fueled by something other than cooling of the hot gas, such as gas-rich mergers (Barnes & Hernquist 1991; Hopkins et al. 2006), it would be reasonable to think that the trends of AGN- and starburst-hosting BCG fractions would be different from the trend of the CC fraction. If this observed increase in AGN activity is linked to gas-rich mergers rather than ICM cooling, we would expect to see an increase in scatter in the P_cav versus L_cool relation (Rafferty et al. 2006) at z > 1.

Another possible scenario to explain the trend of high AGN-hosting BCG fraction at high redshift has to do with cluster mergers. It has been shown both in simulations (Fakhouri et al. 2010) and observations (McDonald et al. 2017) that the cluster merger rate is significantly higher at high redshift. Major mergers between two clusters have the potential of disrupting a tightly regulated AGN feedback loop and promote black hole accretion and star formation by potentially increasing the local turbulence of the system. This is consistent with the turbulent picture in the precipitation model for AGN feedback, called "chaotic cold accretion (CCA)," which states that turbulence is a key component in driving nonlinear thermal instability and extended condensation (Voit et al. 2015; Gaspari et al. 2020). Turbulent forcing can help stimulate precipitation and condensation by increasing the velocity dispersion of the ambient medium, resulting in more black hole accretion and star formation (Voit 2018). With a more energetic environment in the early universe, it is reasonable to assume that the turbulence will be higher at high redshift, resulting in higher black hole accretion rates. The recent discovery of CHIPS 1911 + 4455, a merging galaxy cluster with a massive starburst at the center, provides strong evidence that mergers can indeed increase star formation (Somboonpanyakul et al. 2021a, 2021b). With the development of next-generation X-ray observatories such as Athena and Lynx, we will be able to directly measure motions in the hot gas and determine whether mergers of groups/clusters can boost cooling via an increase in turbulence.

Last, Figure 5 shows that the BCG of the Phoenix cluster remains the most extreme AGN in the entire SPT-SZ sample, which covers a 2500 deg² area and spans all redshifts. In combination with the recent work from the CHiPS survey (Somboonpanyakul et al. 2021a), which confirmed that the Phoenix cluster hosts the most extreme BCG with the strongest CC at z < 0.7, the runaway cooling phase, as we have seen in the Phoenix cluster, is indeed extremely rare.

In this work, we present results for the mid-IR colors of BCGs in SPT-selected galaxy clusters at 0 < z < 1.3. This study allows us to track the evolution of BCG properties over ∼9 Gyr of cluster growth. In particular, we focus our work on black hole accretion in BCGs, which turns these central galaxies into bright AGNs. Our findings are summarized as follows:

• 1.  
  Assuming a single-burst stellar population model, we find statistically significant evidence ( > 99.9%) for a mid-IR excess in high-redshift BCGs compared to low-redshift BCGs, suggesting an increase with redshift in the fraction of AGN-hosting BCGs in galaxy clusters over 0 < z < 1.3. For the lower-redshift bins (z < 0.6), an increase is not statistically significant, and the results are compatible with the noise level. On the other hand, we see an increase in the fraction of BCGs with AGNs in high-redshift bins (z > 0.6), similar to what others have found in previous works (Hlavacek-Larrondo et al. 2013; McDonald et al. 2016; Bîrzan et al. 2017).
• 2.  
  We show that our results are consistent with both the evolution of the fraction of AGNs in field galaxies (Silverman et al. 2009; Haggard et al. 2010) and the fraction of starburst BCGs (Webb et al. 2015; McDonald et al. 2016), suggesting that the reason for the evolution of both AGN and starburst fraction could come from the fact that more cold gas is available in the early universe. This should lead to a higher level of gas-rich mergers in BCGs, which could fuel both AGN activity and star formation at the center of clusters. Some caveats about the direct comparison between cluster and field galaxies remain and range from selection criteria to the evolution of AGN luminosity.
• 3.  
  Another possible explanation for the increase in the fraction of AGN-hosting BCGs with redshift could be a higher level of local turbulence from dynamically active galaxy clusters at high redshift, leading to elevated cooling and subsequent black hole accretion. However, for this scenario, it is difficult to explain the similarity to the trends in field galaxies.
• 4.  
  We do not see any additional cluster with a BCG that is as extreme in the mid-IR color as the Phoenix cluster. In other words, the Phoenix cluster likely hosts the most extreme central AGN in the SPT sample.

An enhancement of AGN activity in BCGs at high redshift compared to low redshift, which is similar in magnitude to the increase observed in field galaxies and cluster members, suggests that this increased AGN activity is not related to cooling flows, but rather to accretion of gas-rich satellites at early times. These accretion events ought to throw off the precise cooling/feedback balance in the centers of clusters that responsible for preventing runaway cooling flows, which leads to a less tightly regulated feedback loop at early times. Further studies with deeper and higher angular resolution mid-IR imaging, such as the upcoming James Webb Space Telescope (JWST; Gardner et al. 2006), will be required to better understand the evolution of AGN feedback and its impact on galaxy clusters.

T.S. and M.M. acknowledge support from the Kavli Research Investment Fund at MIT, and from NASA through Chandra grant GO5-16143.

F.R. acknowledges financial support provided by NASA through SAO Award Number SV2-82023 issued by the Chandra X-Ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060.

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

The South Pole Telescope is supported by the National Science Foundation through grant PLR-1248097. Partial support is also provided by the NSF Physics Frontier Center grant PHY-1125897 to the Kavli Institute of Cosmological Physics at the University of Chicago, the Kavli Foundation, and the Gordon and Betty Moore Foundation grant GBMF 947.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientìfico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Enérgeticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Fìsica d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians-Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the OzDES Membership Consortium, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University.

Facilities: Pan-STARRS - , Chandra X-ray Observatory (ACIS) - , Hubble Space Telescope (ACS/WFC) - , the Nordic Optical Telescope (ALFOSC). -

Software: astropy (Astropy Collaboration et al. 2018), CIAO (Fruscione et al. 2006), pandas (McKinney 2010), seaborn (Waskom et al. 2016).