Cosmic evolution of the incidence of Active Galactic Nuclei in massive clusters: Simulations versus observations

It is well established that the main sites of galaxy formation in the Universe are halos of dark matter. These provide the necessary gravitational potential for the various baryonic physical processes to act and form the luminous structures (i.e. galaxies) we observe. Among the different methods proposed in the literature for associating galaxies (i.e. luminous baryonic matter) with dark matter halos, the semi-empirical approach of abundance matching offers a number of advantages. With relatively small number of parameters, this approach can successfully reproduce observed properties of galaxies such as their stellar masses or SFRs. In the basic implementation of abundance matching it is assumed that most massive halos are associated with the most massive galaxies. This approach yields a relation between dark matter mass and stellar mass as a function of redshift that is in reasonable agreement with observational results (e.g. occupation number, two point correlation function or cross bias, see Vale & Ostriker, 2004; Kravtsov et al., 2004; Tasitsiomi et al., 2004). Recent implementations of abundance matching techniques include an increasing level of complexity in the way halos are associated with baryonic mass and galaxies. For example, the halo mass vs. stellar mass relation is parameterized by analytic functions allowing for intrinsic scatter (e.g. Behroozi et al., 2010; Moster et al., 2010), baryonic process such as star-formation in galaxies are modeled using information on the accretion/merger history of halos, diverse observational results (e.g. stellar mass functions, galaxy clustering properties) are used to tune the various model parameters and produce realistic mock galaxy catalogues out to high redshift (Behroozi et al., 2019; Moster et al., 2018).

In this work we use the UniverseMachine data release 1²²2https://www.peterbehroozi.com/data.html (Behroozi et al., 2019) implemented for the MultiDark PLanck2 (MDPL2, Klypin et al., 2016) dark matter N-body simulation. We choose to use the MDPL2 because it is one of the largest volume, high resolution and public cosmological simulations. It has a box size of 1 000 Mpc/$h$, a mass resolution of $1.5\times 10^{9}\,M_{\odot}/h$ and 3 840³ ($\sim$57$\cdot$10⁹) particles. Individual dark matter halos in the MDPL2 are identified using Rockstar (Behroozi et al., 2013a). This is a state-of-the-art halo finder that uses both the 6-dimensional phase-space distribution of dark matter particles and temporal information to identify bound structures, i.e. dark matter halos. Rockstar is efficient in detecting and measuring the properties of both the largest collapsed structures (parent haloes) and sub-structures within them (satellites haloes). The evolution of haloes through cosmic time is tracked in the form of merger trees computed by the code Consistent Trees (Behroozi et al., 2013b). In this work we consider only dark matter halos with at least 100 times the MDPL2 mass resolution, i.e. $M_{\rm peak}>1.5\times 10^{11}\,M_{\odot}/h$. This limit ensures that the inferred properties of dark matter halos, such as their position and total mass are not affected by the finite resolution of the simulations.

UniverseMachine assigns stellar masses (and hence galaxies) to dark matter halos by parameterizing the star formation history (SFH) of individual halos. The SFR in a halo is assumed to be a function of the depth of the halo’s potential well, its assembly history and cosmic time. The maximum circular velocity, $v_{\rm max}$, is used as a proxy of the depth of the potential well. The $v_{\rm max}$, corresponds to the circular velocity of the halo when it reaches its historical maximum mass ($M_{\rm peak}$ parameter in the MDPL2 catalogues). The halo assembly history is parametrized by the $v_{\rm max}$ variations ($\Delta v_{\rm max}$) across cosmic time. UniverseMachine therefore assumes a parametric analytic function SFR$(v_{\rm max},\Delta v_{\rm max},z)$ to determine the SFR for each halo across cosmic time. Integrating the SFR along the assembly and merger history of a galaxy it is then possible to determine the corresponding stellar mass. The parameter space of SFR$(v_{\rm max},\Delta v_{\rm max},z)$ function is explored in an iterative manner by estimating at each step observables (stellar mass functions, UV luminosity functions, the UV–stellar mass relation, specific and cosmic SFRs, galaxy quenched fractions, galaxy auto-correlation functions and the quenched fraction of central galaxies as a function of environmental density) and comparing them with observations at different redshfits. A Monte Carlo Markov Chain (MCMC) approach is used to sample the model parameter space and yield posteriors for the model parameters.

The end product of UniverseMachine are catalogues of dark matter halos, each of which is assigned a galaxy stellar mass and a SFR. By construction the galaxy population is consistent with observations, including the stellar mass function at different redshift, the evolution of the SFR density of the Universe and the Main Sequence of star-formation. In the following we use the “observed” UniverseMachine values for the stellar mass and SFR of mock galaxies. These are estimated by adding systematic errors (also free parameters in UniverseMachine) to the “true” values to account for observational effects (e.g. Eddington bias). We note however, that our final results and conclusions are not sensitive to this choice. For dark matter haloes we use virial values as defined by Bryan & Norman (1998) for mass and radius.

The assignment of AGN to the UniverseMachine galaxies is also based on empirical relations that associate the probability of a supermassive black hole accretion event to the properties of its host. The relevant observable is the specific accretion rate, $\lambda_{\rm{sBHAR}}\propto L_{X}({\rm 2-10\,keV})/M_{\star}$. In this definition $L_{X}(\rm 2-10\,keV)$ is the AGN X-ray luminosity in the $2-10$ keV band and $M_{\star}$ is the stellar mass of the parent galaxy. The specific accretion rate provides an estimate of how much X-ray luminosity is emitted by the AGN per unit stellar mass of the host galaxy. In this work we choose to scale the ratio $L_{X}({\rm 2-10\,keV})/M_{\star}$ as

Large multi-wavelength observational programs have enabled the estimation of stellar masses, X-ray luminosities and hence $\lambda_{\rm{sBHAR}}$ for large samples of AGN (Aird et al., 2012; Bongiorno et al., 2012; Schulze et al., 2015; Georgakakis et al., 2017). These observations made possible the determination of the fraction of galaxies at fixed stellar mass that host an accretion event with specific accretion rate $\lambda_{\rm{sBHAR}}$. These fractions can then be turned into specific accretion rate probability distribution functions, P($\lambda_{\rm{sBHAR}}$), which describe the probability of an accretion event with parameter $\lambda_{\rm{sBHAR}}$ in a galaxy. Recent observational studies have measured the specific accretion rate distribution as a function of redshift and host galaxy properties such as stellar mass and SFR (Aird et al., 2018; Georgakakis et al., 2017). In this work we use these two independent estimates of the specific accretion rate distribution.

Georgakakis et al. (2017) combined a number of extragalactic X-ray survey fields with multi-wavelength data to construct a non-parametric model of the specific accretion rate distribution. Their methodology required that the convolution of the P($\lambda_{\rm{sBHAR}}$) with the galaxy stellar mass function yields the observed number of X-ray AGN in bins of luminosity, redshift and stellar mass. Aird et al. (2018) started with a sample of near-infrared selected galaxies for which stellar masses and SFRs were estimated. X-ray observations were then used to extract the X-ray photons at the positions of individual galaxies. These were then fed into a flexible Bayesian mixture model to determine in a non-parametric manner the corresponding specific accretion rate distribution of star-forming and quiescent galaxies as a function of stellar mass and redshift. Despite differences in the methodology the Georgakakis et al. (2017) and Aird et al. (2018) constraints on the specific accretion rate distribution are in good agreement (see Georgakakis et al., 2017). Both Georgakakis et al. (2017) and Aird et al. (2018) measured P$(\lambda_{\rm{sBHAR}})$ as a function of redshift out to $z\approx 3$. Figure 1 graphically shows examples of the specific accretion rate distributions used in our analysis.

Using these distributions we associate AGN X-ray luminosities to galaxies in the UniverseMachine catalogues. This process is done in a probabilistic approach. For each mock galaxy with stellar mass $M_{\star}$ and redshift $z$ the corresponding specific accretion rate distributions from either Georgakakis et al. (2017) or Aird et al. (2018) are sampled to draw random $\lambda_{\rm{sBHAR}}$. These are then used to assign X-ray luminosities to individual galaxies by inverting Equation (1). Aird et al. (2018) provides separate P$(\lambda_{\rm{sBHAR}})$ for star-forming and passive galaxies. In this case we split the UniverseMachine galaxies into these two classes using the relation of Aird et al. (2018)

The extragalactic survey fields used by Georgakakis et al. (2017) and Aird et al. (2018) are dominated by low density regions of the cosmic web. Groups and cluster of galaxies, although present in these samples, are subdominant simply because of the form of the halo mass function and the relevantly small FOV of most of the survey fields used. The derived specific accretion rate distributions are therefore representative of the field galaxy population, i.e. those outside massive groups or clusters of galaxies. For this reason in what follows we refer to the predictions of the model as “field expectation”. The adopted specific accretion rate distributions are agnostic to the parent halo mass of individual galaxies, hence the the zero-order assumption of the model that the incidence of AGN in galaxies is independent of environment.

The final products of the AGN seeding process are MDPL2 cosmological boxes at fixed redshift with galaxies (from UniverseMachine) and AGNs (from random sampling of the $\lambda_{\rm{sBHAR}}$ distribution). Figure 2 graphically demonstrates that our semi-empirical methodology by construction reproduces the halo mass function of dark matter haloes, the stellar mass function of galaxies and the AGN X-ray luminosity function. It is also demonstrated that this approach produces AGN mocks with the large-scale clustering ($\ga 1$ Mpc) consistent with observations (Georgakakis et al., 2019; Aird & Coil, 2021). Moreover, the AGN duty cycle, defined as the probability of galaxies above a given stellar mass hosting an AGN above a given accretion luminosity, are inherent in the derivation of specific accretion rate distributions above. As a result our AGN and galaxy mocks are consistent with independently derived determinations of the AGN duty cycles (e.g. Goulding et al., 2014). Put differently, the stellar mass function of the mock AGN host galaxies at fixed X-ray luminosity threshold is consistent with observational constraints (Georgakakis et al., 2011, 2017).

It is noted that the above methodology for seeding galaxies and halos with AGN is similar to that proposed by Shankar et al. (2020) and Allevato et al. (2021) for generating mock AGN samples based on the semi-empirical approach. In these studies satellites and central galaxies/AGN can be treated separately by changing their relative duty cycles. In our approach there is no distinction between the two, i.e. it is assumed that both central and satellites are described by the same duty cycle.

The comparison of the predictions from the simulations with the observations is following the principles of forward modeling. For that the UniverseMachine boxes need to be first projected onto the sky plane to mimic real observations. We assume a box with a XYZ Cartesian coordinate system. This is offset along the Z-axis by the co-moving distance $D_{c}(z)$ at redshift $z$. The observer is placed at $\rm{Z}=0$ Mpc$/h$ and on the XY plane. The line-of-sight angle of the observer relative to every galaxy in the box is then estimated. This angle can be split into a right ascension and declination on the unit sphere. The redshift of each object corresponds to its co-moving distance from the observer. The finite FOV of real observations can also be imposed by defining a sight-line from the observer to a given light-cone direction, estimating the angular distances of all mock galaxies relative to this direction and then rejecting the ones with angular distances larger than the adopted FOV.

For the analysis presented in this paper we construct light-cones in the vicinity of clusters. We consider three redshifts $z=0.2,0.75$ and 1.25, which correspond to the mean redshifts of the cluster samples presented by Martini et al. (2009); Martini et al. (2013). We select UniverseMachine boxes with scale factor³³3Defined as $a=1/(1+z)$. 0.4505, 0.5747 and 0.8376 that approximately correspond to each of the redshifts above. For a given box a massive dark matter halo is selected (see Section 3.4 for details) and is placed at a co-moving distance $D_{c}(z_{c})$ from the observer, where $z_{c}$ is the box redshift, i.e. one of 0.2, 0.75 and 1.25. The light-cone to the cluster is then constructed with an opening angle that is defined by the user (see next section for details). The end product are dark matter haloes, mock galaxies and AGN projected on the sky that mimic real observations.

This section describes how observational selection effects are implemented into our simulations to allow comparison with the results of Martini et al. (2009); Martini et al. (2013) on the fraction of AGN in galaxy clusters. The characteristics of the observations we attempt to mimic can be grouped into three broad categories that relate to the richness/mass of the cluster sample, the galaxy/AGN cluster membership criteria and the apparent brightness or stellar mass of the galaxy/AGN sample. Below we discuss each of them in detail.

We first explore the sensitivity of the estimate AGN fractions to observational selection effects. Figure 4 plots the AGN duty cycle of our semi-empirical model as a function of host galaxy stellar mass limit. The former quantity is defined as the fraction of mock AGN above a given X-ray luminosity threshold among galaxies more massive than a given stellar mass limit (x-axis of Figure 4). In this exercise all galaxies in a given UniverseMachine box are used independent of halo mass. We iterate that the AGN fractions plotted in Figure 4 are consistent with independent observations that use the stellar mass function of galaxies and AGN hosts to infer duty cycles (see Georgakakis et al., 2017). Figure 4 shows that the AGN duty cycle is sensitive to both the X-ray and the host galaxy stellar mass thresholds. There is a general trend of increasing AGN fraction with decreasing X-ray luminosity. This is because less luminous AGN are more common than high accretion luminosity events (i.e. the X-ray luminosity function). Also the AGN duty cycle at fixed luminosity drops with decreasing stellar mass below about $10^{11}\,M_{\odot}$. This is because in our implementation lower stellar mass hosts require higher $\lambda_{\rm{sBHAR}}$ to produce AGN above a given luminosity cut. However, higher specific accretion rates are less likely. This is demonstrated by the form of the specific accretion rate distribution plotted in Figure 1, which strongly decreases with increasing $\lambda_{\rm{sBHAR}}$. For stellar masses $\ga 10^{11}\,M_{\odot}$ the duty cycle curves of Figure 4 either increase, decrease or remain nearly flat with increasing stellar mass. This is related to the stellar mass dependence of the specific accretion-rate distributions of Georgakakis et al. (2017) and Aird et al. (2018) used in our analysis. It is also worth noting that the differences between the two specific accretion rate distributions models above are stronger for the lowest redshift panel ($z=0.2$). This is related to the small sample of low redshift AGN in the Georgakakis et al. (2017) and Aird et al. (2018) studies. In any case, the important point of Figure 4 is that AGN fractions are sensitive to the choice of the stellar mass and X-ray luminosity thresholds, i.e. the observational selections effects. This emphasises the importance of carefully treating this issue, as described in Sections 3.4.3.

For completeness we also show in Figure 5 how the fraction of AGN varies with parent halo mass. UniverseMachine galaxies (central or satellites) above a given stellar mass threshold are grouped by parent halo mass. At a fixed halo mass the fraction of galaxies that host AGN above a given accretion luminosity limit is estimated and plotted. As expected the resulting curves are nearly flat with halo mass since the explicit assumption of the semi-empirical model construction is that the probability of a galaxy hosting an accretion event is agnostic to halo mass. Nevertheless, the strong correlation between halo and stellar mass as well as the X-ray luminosity cut imprint systematic trends onto the curves of Figure 5. These are manifested for example, by the increasing AGN fraction toward lower halo masses. This is more pronounced for the curves with a high stellar mass cut, $\log M_{\star}/M_{\odot}>10.8$. This is because this threshold essentially removes a large number of lower mass galaxies, which are typically found in low mass halos. Additionally the form of the specific accretion rate distributions dictates that more massive galaxies are more likely to host AGN above a fixed accretion luminosity threshold. The net effect is the observed increase in the AGN fraction toward the low halo mass end in the case of the higher stellar mass threshold in Figure 5.

Having demonstrated the strong impact of selection effects on the calculation of AGN fractions, we next turn to the comparison of the our model predictions with the observed fractions of AGN in massive galaxy clusters. Figure 6 plots the fraction of AGN among cluster member galaxies as a function of redshift. The observational results of Martini et al. (2009); Martini et al. (2013) at mean redshifts $z=0.2$, 0.75 and 1.25 are compared with the predictions of our semi-empirical model using either the Aird et al. (2018) or the Georgakakis et al. (2017) specific accretion rate distributions. The model predictions for the mass- and magnitude-limited samples are presented in different panels. Cluster AGN fractions are estimated for two luminosity thresholds, $L_{X}(\rm 2-10\,keV)>10^{43}\,\rm{erg\,s^{-1}}$ and $>10^{44}\,\rm{erg\,s^{-1}}$ indicated by different colors. The uncertainties assigned to the model fractions are determined using bootstrap resampling. For each cluster a total of 10 AGN realisations are generated by repeating the seeding process of Section 3.2 10 times (re-seeded samples). For a given cluster these 10 realisations differ in their AGN populations because of the stochastic nature of the seeding process. This results in an extended cluster sample that is 10 times larger than the original, i.e. 3880, 1570, 180 for $z=0.2$, 0.75 and 1.25 respectively. At fixed redshift clusters are drawn with replacement from the extended sample to generate a total of 100 sub-samples which are used to determine the 68% confidence interval around the mean. These confidence intervals are the errobars of the model predictions plotted in Figure 6. They represent the uncertainty of the mean expected fraction of AGN per cluster.

The different model flavors broadly yield consistent results at fixed redshift and X-ray luminosity threshold. At the lowest redshift bin however, differences are apparent between the AGN fractions at $L_{X}(\rm 2-10\,keV)>10^{44}\,erg\,s^{-1}$ for the models using the Aird et al. (2018) and Georgakakis et al. (2017) specific accretion rate distributions.

These differences are ultimately linked to the observationally measured specific accretion rate distributions and their dependence on the stellar mass of the AGN hosts. The probability of high mass galaxies hosting an accretion event is lower in the Aird et al. (2018) specific accretion rate distributions compared to the Georgakakis et al. (2017) ones. As a result the AGN fractions predicted by the model flavor that uses the Georgakakis et al. (2017) specific accretion rate distributions are higher. This is primarily because of luminous AGN assigned to the central galaxies of the clusters. About 42 out 388 ($\sim$11%) of the mock clusters have a central galaxy with assigned AGN luminosity $L_{X}(\rm 2-10\,keV)>10^{44}\,\rm{erg\,s^{-1}}$. Such a high incidence AGN is inconsistent with observational constraints on the X-ray properties of Brightest Cluster Galaxies in local massive clusters (Yang et al., 2018). This is a limitation of the Georgakakis et al. (2017) observationally determined specific accretion rate distributions.

The fractions predicted by the semi-empirical model flavor with the apparent magnitude selection effects (right panel in Figure 6) does not include the contribution of AGN emission to the mock galaxy SED. In the case of unobscured AGN this contribution may dominate over the host galaxy stellar component in the optical/mid-infrared bands, particularly at high accretion luminosities ($L_{X}\rm\ga 10^{44}\,erg\,s^{-1}$). The estimated model fractions may therefore be underestimated, because a higher fraction of mock AGN would be brighter than the adopted apparent magnitude cut, if their contribution to the mock galaxy SED was modeled. We estimate nevertheless that this effect is small. We modify the methodology of Section 3.4 by including all AGN cluster members above the luminosity limits $L_{X}(\rm 2-10\,keV)>10^{43}\,\rm{erg\,s^{-1}}$ or $>10^{44}\,\rm{erg\,s^{-1}}$ in the calculation of the corresponding AGN cluster fractions, irrespective of the apparent magnitude of their host galaxy. This is equivalent to assuming that all of the mock AGN are unobscured and therefore their apparent magnitudes dominate that of their hosts. This approach nevertheless increases the estimated fraction at any redshift by $\la 0.1$ dex and therefore does not impact our results and conclusions.

A striking result in Figure 6 is the higher observed fraction of AGN in massive clusters at $z=1.25$ compared to the model predictions. The largest discrepancy of nearly 1 dex is for powerful AGN with $L_{X}(\rm 2-10\,keV)>10^{44}\,erg\,s^{-1}$. An enhanced fraction is also observed for moderate luminosity AGN, $L_{X}(\rm 2-10\,keV)>10^{43}\,erg\,s^{-1}$, but the fact that the observations can only place a lower limit does not allow firm conclusions on the amplitude of the effect. These findings can be attributed to systematic differences between field and cluster either in the AGN specific accretion rate distributions or the stellar mass function of the galaxy population. The former option could be interpreted as evidence for environmental dependence of the AGN triggering efficiency. The latter would indicate differences in the galaxy populations as a function of position on the cosmic web.

Observations indicate that the shape of the total (i.e. independent of galaxy type or SFR) stellar mass function of galaxies does not strongly depend on environment out to $z\approx 1.5$ (e.g. Vulcani et al., 2013; Nantais et al., 2016). This behaviour is also reproduced in the UniverseMachine semi-empirical model. This is demonstrated in Figure 7 that compares the (average) stellar mass function of mock galaxies in the same massive clusters used in our analysis with that of all galaxies in the UniverseMachine box. The mass functions in Figure 7 are normalised so that their integral yields the same stellar mass density. This allows direct comparison of the mass function shapes, which are remarkably similar between massive clusters and the full box. The latter is dominated by field galaxies (i.e. not associated with massive halos) and by construction reproduces the observed stellar mass function at different redshifts estimated using extragalactic survey fields (e.g. Muzzin et al., 2013; Ilbert et al., 2013; Moustakas et al., 2013). The construction of AGN mocks using the seeding process described in section 3.2 is primarily sensitive to the shape of the galaxy mass function. The evidence above therefore suggests that the difference between observations and semi-empirical model predictions in Figure 6 cannot be attributed to systematic variations of the total stellar mass function with environment.

We acknowledge differences between field and cluster mass functions for star-forming and quiescent galaxies (e.g. Vulcani et al., 2013; Nantais et al., 2016; Papovich et al., 2018; van der Burg et al., 2020), in the sense that dense environments host a larger fraction of quenched galaxies. Nevertheless, such variations are second order effect in our empirical AGN-seeding model. This is shown in Figure 6, where the predictions using the Aird et al. (2018) specific accretion rate distribution model that includes star-formation dependence and those based on the Georgakakis et al. (2017) specific accretion rate distribution (no SFR-dependence) are similar at $z=1.25$.

The higher fraction of AGN in $z\ga 1$ clusters compared to the field expectation in Figure 6 contradicts the results of Martini et al. (2013), who report similar fractions. The field AGN fraction of Martini et al. (2013) is estimated from annular regions centered on the clusters with inner and outer radii of 2 and 6 arcmin respectively. These regions may include filaments, infalling groups and generally dense structures associated with the nodes of the cosmic web where the cluster is found. They may therefore not be entirely representative of the true field. Relevant to this point are recent results by Koulouridis & Bartalucci (2019) who studied the radial distribution of massive galaxy clusters at $z\approx 1$. They report a statistically significant overdensity of X-ray selected AGN within the infall cluster region at a project distance of $\sim 2-2.5\,R_{500}$ relative to the cluster center. This interval lies within the region used by Martini et al. (2013) to determine their field AGN fractions. The enhanced number of AGN found by Koulouridis & Bartalucci (2019) could bias high the field AGN fractions estimated by Martini et al. (2013). Our findings are consistent with the higher fraction of infrared selected AGN in massive clusters of galaxies at $z\ga 1$ reported by Alberts et al. (2016). We caution nevertheless that the selection function of that sample is very different from that of Martini et al. (2013). Alberts et al. (2016) identifies AGN by fitting model templates to the multi-wavelength SEDs of mid-infrared selected cluster members.

It is interesting to speculate on specific accretion rate distributions that reproduce the Martini et al. (2013) AGN fractions in massive clusters of galaxies at $z\approx 1.25$ shown in Figure 6. A distribution is required that produces more luminous AGN compared to the current model. There are clearly many functional forms that could achieve that. For simplicity we approach this problem by assuming a Gaussian for the specific accretion rate distribution with parameters (mean, scatter) free to vary. We caution that this problem has a broad range of non-unique solutions. This is ultimately related to the limited observational constraints that are not sufficient to break degeneracies among model parameters. There are essentially only two data points, one of which an upper limit, on the cluster AGN fraction at $z\approx 1.25$ shown in Figure 6. It is nevertheless, instructive to explore parameter combinations (mean, scatter) that yield AGN fractions consistent with the observations. In practice we assume a Gaussian specific accretion rate distribution with a given mean/scatter that is independent of stellar mass or SFR. This is applied to UniverseMachine galaxies (i.e. similar to Section 3.2) to produce light-cones of clusters (Section 3.3) on which selection effects are applied (Section 3.4). The resulting AGN factions are then required to be within the $1\,\sigma$ uncertainty of the Martini et al. (2013) $L_{X}>\rm 10^{44}\,erg\,s^{-1}$ data point or larger than the lower limit for $L_{X}>\rm 10^{43}\,erg\,s^{-1}$ AGN in Figure 6. The general trend is that broader distributions (larger scatter) require lower means $\lambda_{sBHAR}$ to reproduce the data. We also find that Gaussians with mean of $\lambda_{sBHAR}\ga 10^{-1}$ produce too many luminous AGN for any scatter value and are therefore not allowed by the observations. Similarly a mean value of $\lambda_{sBHAR}\la 10^{-3}$ produce too few luminous AGN for any scatter value and are also rejected.

Figure 8 (left panel) shows a selected number of Gaussian specific accretion rate distributions which produce AGN fractions consistent with the observations. These distributions have more power at intermediate specific accretion rates ($10^{-2}\la\lambda_{sBHAR}\la 10^{-1}$) compared to the Georgakakis et al. (2017) and Aird et al. (2018) models that represent the field AGN specific accretion rate distributions. The Gaussian specific accretion rate models plotted in Figure 8 make different predictions on the number of cluster AGN as a function of accretion luminosity. This is demonstrated in the middle panel of Figure 8, which plots the predicted cluster AGN X-ray luminosity functions for the different Gaussian specific accretion rate distribution models. The cluster XLF prediction is different in both shape and normalisation from the field one, also plotted in Figure 8. Future observations that provide a broad luminosity baseline and sufficient AGN number statistics can help constrain the models by e.g. directly measuring the XLF as a function of environment. We have also confirmed that even if massive clusters in the UniverseMachine box are assigned the higher normalisation XLFs shown in Figure 8, the total AGN XLF averaged across environments is consistent with observations. This is because massive clusters are rare as a result of the shape of the Halo Mass Function and therefore have a minor contribution to the mean XLF of the UniverseMachine box.

Finally, the difference between field and cluster specific accretion rates plotted in the left panel of Figure 8 can be the result of varying black hole Eddington ratio distributions and/or duty cycles as a function of environment. These are the physical quantities that convolve to yield the specific accretion rate distribution (e.g. Shankar et al., 2020; Allevato et al., 2021). The analysis presented in this paper cannot identify which of the two quantities (Eddington ratio or duty cycle) is primary driving the differences in the AGN fraction between field and cluster environments. Further work is needed to address this issue and associate the observed differences to physical quantities directly related to the accretion flow onto the supermassive black hole.

Contrary to the results at $z=1.25$ discussed above, Figure 6 shows that at lower redshift, $z=0.75$, the fraction of AGN in massive clusters is consistent with the model predictions (i.e. field expectation, see Section 3.2) within the observational data uncertainties. This conclusion does not strongly depend on the details of the semi-empirical modelling, e.g. which specific accretion rate distribution is adopted for seeding galaxies with AGN or the type of observational selection effects (mass vs. apparent magnitude cut) applied to the mock sample. At even lower redshift, $z=0.2$ our analysis tentatively suggests a paucity of AGN in cluster galaxies compared to the field. This is mainly driven by the higher (factor 2–3) fraction of AGN with $L_{X}(\rm 2-10\,keV)>10^{43}\,erg\,s^{-1}$ predicted by the model compared to the observations. For more luminous AGN, $L_{X}(\rm 2-10\,keV)>10^{44}\,erg\,s^{-1}$, no firm conclusions can be made because the observations only provide an upper limit to the AGN fraction in clusters. We also caution that at this luminosity cut the semi-empirical model that uses the Georgakakis et al. (2017) specific accretion rate distribution is biased high. In this model flavor a large fraction ($\approx 10$%; see Section 3.4) of the massive central galaxies is assigned powerful AGN. This fraction is much higher than the observed incidence of luminous AGN ($\approx 1-3$%) among the Brightest Cluster Galaxies Yang et al. (2018). This discrepancy is ultimately related to the stellar mass dependence of the Georgakakis et al. (2017) empirical specific accretion rates. A stronger such dependence exists in the Aird et al. (2018) specific accretion rate distribution. As a result this model predicts much lower $L_{X}(\rm 2-10\,keV)>10^{44}\,erg\,s^{-1}$ AGN fractions in $z=0.2$ clusters.

Overall the evidence above points to a differential redshift evolution of the incidence of AGN in clusters relative to the field. Massive structures at $z\ga 1$ are found to be more efficient in triggering accretion events onto supermassive black holes compared to less dense regions. Such an environmental dependence is not present at $z\approx 0.75$ and possibly inverses at low redshift, $z\approx 0.2$, in the sense that clusters likely become less active regions for black hole growth compared to the field.

Eastman et al. (2007) also find evidence for differential evolution of the AGN population as a function of environment. They compared the fraction of X-ray selected AGN in cluster of galaxies between $z\approx 0.2$ and 0.6. Their analysis suggests an accelerated evolution of the AGN population in dense environments compared to the field between these redshifts. X-ray observations of individual proto-clusters at higher redshift $z\ga 2$ find that AGN are more common among galaxies in these environments compare to the field (Lehmer et al., 2009; Digby-North et al., 2010; Krishnan et al., 2017). These trends suggests that the probability of a galaxy hosting an accretion event depends on its small-scale environment, $\rm<1\,Mpc$, in agreement with our findings.

It is also interesting that observational studies of the Halo Occupation Distribution (HOD) of AGN are broadly consistent with this picture. At low-redshift, $z\la 0.2$, the HOD of X-ray AGN is proposed to have satellite fractions that increase with halo mass, albeit less rapidly than the galaxy population (Miyaji et al., 2011; Allevato et al., 2012). This points to a decreasing fraction of AGN relative to galaxies with increasing halo mass, similar to our conclusions for low redshift clusters. At higher redshift however, studies of the quasi-stellar object projected correlation function find that the quasi-stellar object satellite fractions increase with halo mass similar to galaxies (Richardson et al., 2012; Shen et al., 2013). The evidence above is therefore consistent with a picture whereby the fraction of AGN in massive clusters evolves with redshift.

Next we explore different physical mechanisms that could be responsible for the trends discussed in the previous sections. In that respect it is interesting that the redshift evolution of the fraction of AGN in clusters relative to the field bears similarities to the star-formation properties of cluster member galaxies as a function of cosmic time. In the local Universe clusters are dominated by quiescent galaxies and their integrated SFRs are significantly lower than the field expectation (e.g. Chung et al., 2011). This changes however toward higher redshift. The total cluster SFR normalised to halo mass increases with lookback time and possibly catches up with the field at $z>1$ (e.g. Webb et al., 2013; Popesso et al., 2015; Alberts et al., 2016). This is also accompanied by an inversion of the SFR vs. density relation (Butcher & Oemler, 1984), whereby cluster cores at $z\ga 1$ host large fractions of actively star-forming galaxy populations (e.g. Elbaz et al., 2007; Alberts et al., 2016). The trends above could be explained by the increasing fraction of the cold gas content of galaxies toward higher redshift (e.g. Tacconi et al., 2010; Saintonge et al., 2013; Santini et al., 2014; Gobat et al., 2020) and the increasing fraction of galaxy interactions in overdense regions compared to the field. However, it is still unclear if cluster member galaxies at $z\ga 1$ are as gas rich as their field counterparts. Stacking analysis at millimeter wavelengths using ALMA continuum observations find that star-forming galaxies in massive clusters at $z\ga 1$ are on average significantly more deficient in molecular gas relative to the field (Alberts et al., 2022). This is at odds with CO emission-line observations of galaxies in clusters that typically detect molecular gas fractions comparable or higher than the field (Noble et al., 2019; Williams et al., 2022). Despite this discrepancy the general picture emerging from these studies is that star-formation in massive clusters is associated with infalling galaxies during their first passage through the dense region that either manage to retain/replenish their molecular gas (e.g. Kotecha et al., 2022; Vulcani et al., 2018) or have been largely stripped but continue to consume any remaining gas by forming stars before being rapidly quenched (e.g. “delay then rapid” scenario Wetzel et al., 2013) .

It is therefore possible that AGN in clusters are also associated with galaxies that fall for the first time into the deep potential well of the overdensity. Relevant to this point is the discovery of an excess of X-ray AGN at the outskirts (2–3 $R_{500}$) of massive clusters ($M_{500}>10^{14}\,M_{\odot}$) at $z\approx 1$ (Koulouridis & Bartalucci, 2019). This finding points to a direct link between the growth of supermassive black holes and the infall region of massive structures in the cosmic web. Therefore the strong evolution of the AGN fraction in clusters may be linked to the increasing molecular gas content with increasing redshift of the galaxies that fall into massive halos. It is then interesting to speculate on the physical mechanisms that could lead to the differential evolution of the AGN fraction in clusters relative to the field, possible processes include:

1. 1.

   Galaxy interactions. These are expected to be more common in dense environments and the outskirts of clusters. Moreover numerical simulations show that the merging rate of dark matter halos as well as the accretion rate onto them increases with redshift (e.g. Gottlöber et al., 2001; Fakhouri & Ma, 2008; McBride et al., 2009). Therefore the flow of galaxies from the cosmic web onto massive halos is expected to be higher at earlier times therefore leading to enhanced galaxy interaction rates. The Semi-Analytic Model of Gatti et al. (2016) suggests that galaxy interactions dominate relatively luminous AGN $L_{X}\ga 10^{43}\rm\,erg\,s^{-1}$ in massive halos ($\ga 10^{14}\,M_{\odot}$) and low/intermediate redshift, $z\la 1$. Instead, internal processes, such as disk instabilities, become important for lower luminosity AGN in dense environments. At much higher redshift, $z>2$, the relative contribution of galaxy interactions and internal processes in triggering AGN inverses, with disk instabilities dominating in massive halos. This is because of the higher gas and disk fraction of galaxies at these higher redshift (Gatti et al., 2016).

2. 2.

   Ram-pressure. Interaction between the Intra-Cluster Medium (ICM) and the cold gas of the galaxy, can make the latter lose angular momentum (without stripping it), hence facilitating its flow to the nuclear regions where it can accrete onto the supermassive black hole. Simulations have shown that this is possible, if the density of the ICM is not very dense (e.g. in the outskirts of clusters) (e.g. Marshall et al., 2018; Ricarte et al., 2020). Observations of galaxies in local clusters with morphological evidence for ongoing ram-pressure stripping indeed reveal a high incidence ($\sim 27-51\%$) of AGN optical emission-line signatures (Peluso et al., 2022). It is further expected that the ICM is more tenuous toward higher redshift and in the cluster outskirts. This is the regime where the physical conditions are more favorable for the ram-pressure to have a positive impact on AGN triggering. Such a scenario could lead to the observed trend of increasing AGN fractions in clusters relative to the field with increasing redshift. Clusters at $z\approx 0.2$ have dense ICMs and the ram-pressure is sufficiently strong to strip galaxies off their gas even at the outskirts (e.g. Zinger et al., 2018; Arthur et al., 2019) and hence, perhaps suppress powerful AGN. Instead, at higher redshift the conditions may be appropriate for the ram-pressure to have a positive effect and boost the numbers of luminous AGN. The radial distribution of AGN within massive halos could provide further constraints this scenario (e.g. Marshall et al., 2018).

In addition to the processes above, simulations also highlight the potential of cosmic filaments in channeling matter and streams of cold gas deep into the potential well of massive overdensities in the Universe (e.g. Kereš et al., 2005; Dekel & Birnboim, 2006; Kereš et al., 2009), particularly at high redshift, $z\ga 1$. Galaxies that fall into clusters along such filaments are shielded from the hot ICM and can therefore retain at least part of their gas reservoirs even close to the cluster core (Kotecha et al., 2022). This delays the quenching of galaxies by modulating ram-pressure stripping or strangulation and allows the formation of new stars in dense cluster environments. It can be expected that the same process also promotes black hole growth in the galaxies that accrete onto massive halos through filaments. Furthemore Kotecha et al. (2022) argue that this effect is more pronounced at high redshifts $z=1-3$ since at this epoch the cold flow filaments are expected to have a higher temperature contrast relative to the ICM and cluster halos are typically less massive. Such an evolution pattern is consistent with the higher fraction of AGN in massive clusters at $z\ga 1$ compared to lower redshift. Additional processes (e.g. interactions, positive impact of ram-pressure) need to be invoked to explain the difference between the AGN fractions in cluster and field at $z\ga 1$.