Cosmic evolution of the incidence of active galactic nuclei in massive clusters: simulations versus observations

ABSTRACT

1 INTRODUCTION

It is now widely accepted that supermassive black holes are found at the centres of most, if not all, galaxies in the local Universe (Kormendy & Ho 2013). These compact objects are believed to have grown their masses throughout cosmic time via accretion of material from their surroundings (e.g. Soltan 1982; Marconi et al. 2004; Merloni & Heinz 2008). Such accretion events generate large amounts of energy that can be detected as radiation across the electromagnetic spectrum. The astrophysical sources associated with such events are dubbed active galactic nuclei (AGNs; Padovani et al. 2017). Observational campaigns in the last 20 yr aiming at detecting and characterizing large samples of AGNs have painted a comprehensive picture of the cosmological evolution of this population and have provided a quantitative description of the accretion history of the Universe out to high redshift (e.g. Ueda et al. 2014; Aird et al. 2015; Brandt & Alexander 2015). What remains challenging to understand, however, are the physical processes that trigger accretion events on to supermassive black holes and therefore drive the observed black hole growth as a function of redshift. The different supermassive black hole fueling mechanisms proposed in the literature can broadly be grouped into external (ex situ) and internal (in situ) in nature. The latter are related to the secular evolution of galaxies, e.g. disc instabilities (e.g. Hopkins & Hernquist 2006; Gatti et al. 2016), the creation of bars (Cisternas et al. 2015), stellar winds (Ciotti & Ostriker 2007; Kauffmann & Heckman 2009) or the biased collapse of the baryons in the inner region of the halo (e.g. Lapi et al. 2011, 2018), and could lead to gas inflows towards the central regions of the galaxy and feeding of the central black hole. Ex situ processes are those that act on a galaxy from its environment. They include for example, galaxy interactions (Di Matteo, Springel & Hernquist 2005; Gatti et al. 2016), cold gas inflows (Bournaud et al. 2012; DeGraf et al. 2017), or cooling flows in massive clusters (e.g. Fabian 1994). The balance between ex situ and in situ supermassive black hole fuelling processes likely depends, among others, on redshift, position or the cosmic web and intrinsic galaxy properties, such as gas content and structural parameters.

One approach for exploring the relative importance of the diverse mechanisms above in modulating the growth of supermassive black holes is to study the incidence of AGNs in galaxies as a function of e.g. their morphology, star formation rate (SFR) or position on the cosmic web (Brandt & Alexander 2015). Such investigations can shed light on the conditions that promote or suppress accretion events on to the supermassive black holes of galaxies and make inferences on the physics at play. Environmental studies in particular, i.e. how AGNs populate galaxy groups, filaments, clusters, and field, could provide information on the balance between in situ and ex situ process for activating supermassive black holes. This potential have motivated observational studies to characterize AGNs populations in different environments. At low redshift (z ≈ 0.1) there is evidence that the fraction of AGNs in high density regions is lower compared to the field (e.g. Kauffmann et al. 2004; Popesso & Biviano 2006; Koulouridis & Plionis 2010; Lopes, Ribeiro & Rembold 2017; Mishra & Dai 2020). This may indicate the decreasing incidence of mergers in massive clusters (Popesso & Biviano 2006) and/or the impact of processes that strip the gas reservoirs of galaxies and hence, lead to the suppression of their nuclear activity. There are also claims that the AGN radial distribution is skewed to the cluster outskirts relative to the general galaxy population (e.g. de Souza et al. 2016; Lopes et al. 2017). This finding coupled with suggestions that cluster AGN show high velocity dispersion (e.g. Haines et al. 2012; Pimbblet et al. 2013; Lopes et al. 2017) points to a link between accretion events and galaxies that fall on to high density regions from larger scales. Contrary to the findings above there are also observational studies that claim little or no dependence of the AGN fraction on environment at low redshift (e.g. Miller et al. 2003; Haggard et al. 2010; Pimbblet et al. 2013). At least part of the discrepancy is likely related to selection effects. These include the accretion luminosity threshold adopted in the various studies (Kauffmann et al. 2004; Pimbblet et al. 2013), differences between field and cluster environments in the properties of the overall galaxy population (e.g. SFR, morphology) used to determine fractions (e.g. von der Linden et al. 2010; Lopes et al. 2017; Man et al. 2019) or the methods adopted for selecting AGNs (e.g. optical emission lines, X-ray emission, and mid-infrared colours).

Outside the local Universe (⁠|$z\gtrsim 0.1$|⁠) there is evidence that the group/cluster environments become more active in terms of black hole growth. The fraction of AGN in such dense regions increases with increasing redshift (Martini, Sivakoff & Mulchaey 2009; Martini et al. 2013; Bufanda et al. 2017) at a rate that appears to be faster than the field AGN evolution (Eastman et al. 2007). At redshift |$z\gtrsim 1$| the fraction of AGN in clusters is at least as high as the field expectation (Martini et al. 2013; Alberts et al. 2016) suggesting efficient triggering of accretion events. This is possibly associated with the higher incidence of interactions in these environments (e.g. Alberts et al. 2016) and/or the larger cold gas content of galaxies at earlier times (e.g. Tacconi et al. 2010) combined with the impact of the ram pressure experienced by galaxies as they fall into the cluster potential well (Poggianti et al. 2017; Ricarte et al. 2020). The evidence above emphasizes the role of environment for understanding AGN triggering mechanisms and underlines the need to better constrain the redshift at which the cluster AGN fractions are on par with the field or even exceed it (Alberts et al. 2016).

In this work, we revisit the incidence of AGNs in massive clusters of galaxies out to z ≈ 1.25 by developing a semi-empirical modelling approach to interpret observational results from the literature (Martini et al. 2009, 2013). The feature of our modelling methodology is control over systematics and observational selection effects. We use observationally motivated relations to populate massive dark matter haloes extracted from N-body cosmological simulations with AGNs and galaxies. These are then used to mimic observations of clusters of galaxies by including in a realistic manner the relevant selection effects, such as cluster membership definition, flux, or luminosity cuts, etc. These mocks are then compared with real observations to make inferences on the evolution of the AGN fraction in clusters relative to the field expectation. Section 2 presents the observations and selection bias that we attempt to reproduce. Section 3 describes the generation of the mock catalogues and the implementation of the different selection effects into the simulations. The comparison of the semi-empirical model predictions with the observations is presented in Section 4. Finally, Section 5 discusses the results in the context of AGN triggering mechanisms. We adopt a flat ΛCDM cosmology with parameters Ω_m = 0.307, Ω_Λ = 0.693, h = 0.678 consistent with the Planck results (Planck Collaboration XIII 2016).

2 OBSERVATIONS

This work uses the observational measurements of the fraction of X-ray AGN in galaxy clusters presented by Martini et al. (2009, 2013). Typical halo masses of these clusters are few times |$10^{14}\, M_\odot$| (see Section 3.4.1 for more details). In this section, we describe the most salient details of these observations and the corresponding data analysis. Of particular interest to our work are the (i) the definition of cluster membership in the observations and (ii) the magnitude/flux limits that are used to define the galaxy and AGN samples. The inferred AGN fractions strongly depend on these selection effects and it is therefore important to reproduce them in the simulations before comparing with the observations.

Martini et al. (2009) used a sample of 32 massive galaxy clusters out to redshift z = 1.3 with available Chandra X-ray observations. Their low redshift sub-sample consists of 17 clusters at z < 0.4 (mean redshift |$\bar{z}=0.19$|⁠). These 17 clusters include the 10 presented in Martini et al. (2006) and 7 additional ones selected from the Chandra archive to be the nearest most massive clusters with virial radius that fits within the Chandra field of view (FOV). The high redshift sub-sample numbers 15 clusters in the redshift interval z = 0.4–1.3 (average redshift |$\bar{z}=0.72$|⁠). Cluster member candidates are selected within the projected R₂₀₀¹ radius of each cluster in the sample. The number of AGN and galaxy cluster members is determined to the apparent R-band magnitude limit |$m_{R}^{*}(z)+1$|⁠, where |$m_{R}^{*}(z)$| is the break of the R-band luminosity function at the cluster redshift. The latter is estimated assuming that the absolute magnitude break of the luminosity function evolves as |$M_R^{*}(z) = M_R^{*}(z=0)-z$| with |$M_R^{*}(z=0)$| from Christlein & Zabludoff (2003) and early-type galaxy spectral energy distribution for the K-correction. For clusters with a high redshift identification completeness the number of galaxy members is estimated by counting sources with R-band magnitude brighter than |$m_R^{*}(z)+1$| and redshift difference (Δz) relative to the mean cluster redshift (z), |$\Delta z \cdot c \lt 3\, \sigma _v (1+z)$|⁠, where σ_v is the cluster velocity dispersion and c the speed of light. For clusters with limited spectroscopic redshift follow-ups (mostly high-redshift sub-sample) the number of galaxy members is estimated using the cluster richness versus velocity dispersion relation of Becker et al. (2007). This empirical relation is calibrated to yield the number of early-type galaxy cluster members that are more luminous than 0.4L* (i.e. equivalent to |$m_R^{*}(z)+1$|⁠) within the R₂₀₀ radius. AGN cluster members are also selected to have apparent magnitude brighter than |$m_R^{*}(z)+1$| and redshifts that are consistent with |$\Delta z\cdot c \lt 3\, \sigma _v(1+z)$|⁠, i.e. similar to galaxies. The observed number of X-ray AGN cluster members is also corrected for the spectroscopic completeness of each cluster (typically >60 per cent for AGNs). The depth of the Chandra X-ray observations means that AGN samples are complete to hard-band luminosities |$L_X(\rm 2\rm{-}10\, keV) \gt 10^{43}\, erg\, s^{-1}$|⁠. Less luminous X-ray sources suffer incompleteness because of the sensitivity of the Chandra observations and are not used for the estimation of AGN fractions.

The Martini et al. (2013) cluster sample is composed of 13 of the most statistically significant extended X-ray sources detected in the Chandra survey of the Bootes field with spectroscopic identifications in the redshift interval z = 1–1.5 (Eisenhardt et al. 2008). Cluster member candidates, AGN or galaxies, are selected to lie within the projected R₂₀₀ radius and have Spitzer|$3.6\, \mu \rm m$|-band apparent magnitude brighter than |$m_{3.6}^{*}(z)+1$|⁠. The quantity |$m_{3.6}^{*}(z)$| is the break of the |$3.6\, \mu \rm m$| luminosity function at redshift z adopted from Mancone et al. (2010). Both spectroscopic and photometric redshifts are used to determine cluster membership. Sources with spectroscopic redshift within |$\Delta z\cdot c\lt \pm 2\, 000 (1+z)\, \rm km\, s^{-1}$| (i.e. similar to Martini et al. 2009 fixing |$3\sigma _v=2\, 000\, \rm km\, s^{-1}$|⁠) off the cluster redshift are assumed to be members. In the case of photometric redshift estimates this condition is modified so that at least 30 per cent of the photometric redshift probability density function is required to lie within the above redshift interval. All X-ray selected AGN cluster member candidates in the sample of Martini et al. (2013) have spectroscopic redshifts. The AGN sample is complete to the X-ray luminosity |$L_X(\rm 2\rm{-}10\, keV) =10^{44}\, erg \, s^{-1}$|⁠. For less luminous systems, |$L_X(\rm 2\rm{-}10\, keV) \gt 10^{43}\, erg \, s^{-1}$|⁠, Martini et al. (2013) provide lower limits for the AGN cluster fraction.

3 METHODOLOGY

This section describes our approach for generating mock observations of galaxies and AGNs in massive structures of the cosmic web. The starting point of our method is cosmological N-body simulation (e.g. Lemson & Virgo Consortium 2006; Klypin, Trujillo-Gomez & Primack 2011; Klypin et al. 2016) that describe the formation and evolution of dark matter haloes in the Universe under the influence of gravity. These are coupled with empirical relations that associate dark matter haloes with galaxies (galaxy–halo connection). Accretion events associated with supermassive black holes are then painted on top of those galaxies using recent observational results on the incidence of AGN in galaxies (AGN–galaxy connection). The implicit assumption of this latter step is that the probability of galaxies hosting an accretion event does not depend on environment, i.e. halo mass. Light-cones are then generated to mimic real observations of AGNs and galaxies on the cosmic web. These steps above are described in detail in the following sections.

3.1 Galaxy–halo connection: (Sub-)halo abundance matching techniques

It is well established that the main sites of galaxy formation in the Universe are haloes 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 haloes, 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 haloes 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 Kravtsov et al. 2004; Tasitsiomi et al. 2004; Vale & Ostriker 2004). Recent implementations of abundance matching techniques include an increasing level of complexity in the way haloes are associated with baryonic mass and galaxies. For example, the halo mass versus stellar mass relation is parametrized by analytical functions allowing for intrinsic scatter (e.g. Behroozi, Conroy & Wechsler 2010; Moster et al. 2010), baryonic process such as star formation in galaxies are modelled using information on the accretion/merger history of haloes, diverse observational results (e.g. stellar mass functions and galaxy clustering properties) are used to tune the various model parameters and produce realistic mock galaxy catalogues out to high redshift (Moster, Naab & White 2018; Behroozi et al. 2019).

In this work, we use the universemachine data release 1² (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\, \mathrm{M}_\odot \,h^{-1}$| and 3 840³ (∼57 × 10⁹) particles. Individual dark matter haloes in the MDPL2 are identified using rockstar (Behroozi, Wechsler & Wu 2013a). This is a state-of-the-art halo finder that uses both the 6D phase-space distribution of dark matter particles and temporal information to identify bound structures, i.e. dark matter haloes. 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 haloes with at least 100 times the MDPL2 mass resolution, i.e. |$M_{\rm peak}\gt 1.5\times 10^{11}\, \mathrm{M}_\odot \,\mathrm{h}^{-1}$|⁠. This limit ensures that the inferred properties of dark matter haloes, 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 haloes by parametrizing the star formation history (SFH) of individual haloes. 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_max, is used as a proxy of the depth of the potential well. The v_max, corresponds to the circular velocity of the halo when it reaches its historical maximum mass (M_peak parameter in the MDPL2 catalogues). The halo assembly history is parametrized by the v_max variations (Δv_max) across cosmic time. universemachine therefore assumes a parametric analytic function SFR(v_max, Δv_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_max, Δv_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 autocorrelation 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 haloes, each of which is assigned a galaxy stellar mass and an 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.

3.2 AGN–galaxy connection: specific accretion rate distributions

Large multiwavelength observational programs have enabled the estimation of stellar masses, X-ray luminosities and hence λ_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 λ_sBHAR. These fractions can then be turned into specific accretion rate probability distribution functions, P(λ_sBHAR), which describe the probability of an accretion event with parameter λ_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 (Georgakakis et al. 2017; Aird, Coil & Georgakakis 2018). 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 multiwavelength data to construct a non-parametric model of the specific accretion rate distribution. Their methodology required that the convolution of the P(λ_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(λ_sBHAR) as a function of redshift out to z ≈ 3. Fig. 1 graphically shows examples of the specific accretion rate distributions used in our analysis.

Figure 1.

Specific accretion rate distributions that describe the probability of a galaxy hosting an AGN with specific accretion rate λ_sBHAR. The shaded regions correspond to different observational measurements of P(λ_sBHAR). The purple colour shows the Aird et al. (2018) result, where different line styles indicate different galaxy types: star-forming galaxies (solid line) or passive (dashed–dotted). The green colour shows the Georgakakis et al. (2017) constraints on the specific accretion are distribution. All curves correspond to the redshift interval z = 0.5 − 1 and stellar mass interval M_⋆ = 10^10.5–10¹¹ M_⊙. The extent of the shaded regions correspond to the 68 per cent confidence interval around the median (bold curves).

Open in new tabDownload slide

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 zero-order assumption of the model that the incidence of AGNs 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 λ_sBHAR distribution). Fig. 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 (⁠|$\gtrsim 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 AGNs 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).

Figure 2.

Graphical workflow of the semi-empirical modelling to construct mock AGN catalogues based on dark matter N-body simulations. The top panels correspond to a 10 Mpc h⁻¹ slice of a box from MDPL2 cosmological simulation with 1 000 Mpc h⁻¹ side size at the snapshot z = 0.75. The dots represent the positions of dark matter haloes with masses |$M_{\rm peak}\gt 10^{12}\, \mathrm{M}_\odot \,\mathrm{h}^{-1}$| (top left panel), galaxies within these dark matter haloes (top middle panel) and AGN within the same dark matter haloes (top right panel). Only AGN with |$L_X\, (2\rm{-}10\, \rm {keV})\gt 10^{42}\, \rm {erg\, s^{-1}}$| using the specific accretion distribution described in Georgakakis et al. (2017). The construction of AGN mocks proceeds from left to right in this figure. Dark matter haloes (black dots in the top left panel) are populated with galaxies (blue points in the top middle panel), using abundance matching (universemachine; Behroozi et al. 2019). These galaxies are seeded with accretion events following the observationally derived distributions of these events (e.g. Georgakakis et al. 2017; Aird et al. 2018). The key feature of this approach is the reproduction by construction of the predictions and observables shown in the lower set of panels. Lower left panel: halo mass function predicted by theoretical models in orange (HaloMod, https://pypi.org/project/halomod/) and simulations in blue (MDPL2; Klypin et al. 2016). Lower centre panel: stellar mass function where circles represent observations (Ilbert et al. 2013; Muzzin et al. 2013) and the dashed curve the semi-empirical model prediction. Lower right panel: X-ray luminosity function where circles represent observations (Georgakakis et al. 2017) and the curves represent the two independently derived models using the specific accretion rate distributions of Georgakakis et al. (2017, green dotted–dashed) and Aird et al. (2018, purple dashed).

Open in new tabDownload slide

It is noted that the above methodology for seeding galaxies and haloes 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/AGNs 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.

3.3 Light-cones

The comparison of the predictions from the simulations with the observations is following the principles of forward modelling. For that the universemachine boxes need to be first projected on to 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 {\mathit{ 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 sightline 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) and Martini et al. (2013). We select universemachine boxes with scale factor³ 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.

3.4 Selection effects

This section describes how observational selection effects are implemented into our simulations to allow comparison with the results of Martini et al. (2009, 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.

3.4.1 Cluster sample

We define the cluster sample by adopting a minimum virial mass threshold. Martini et al. (2009) provide velocity dispersion for their cluster sample as a measure of their masses. However, the universemachine data set does not include velocity dispersion information and therefore a mapping is required between this parameter and halo mass. For the latter, we adopt the analytical relations presented by Munari et al. (2013) based on N-body simulations. This allows us to associate individual universemachine haloes with a velocity dispersion and then threshold on this quantity to mimic the Martini et al. (2009) cluster sample selection. We choose the lower halo mass limit in such a way that our parent cluster sample reproduces the median velocity dispersion of the Martini et al. (2009) sample. The adopted virial halo mass limits are 5.7 × 10¹⁴ M_⊙ h⁻¹ and 3.6 × 10¹⁴ M_⊙ h⁻¹ for z = 0.2 and 0.75, respectively. For the high-redshift clusters of Martini et al. (2013) there is only scattered information on their halo masses. Literature results suggest masses of few times |$10^{14}\, M_\odot$|⁠, based on dynamical measurements or estimates from X-ray luminosities. For our high-redshift simulations (z = 1.25) we therefore select haloes with virial mass |$\gt 3\times 10^{14}\, \mathrm{M}_\odot\,{\mathrm {h}^{-1}}$| to mimic the cluster sample of Martini et al. (2013). Our results and conclusions are not sensitive to this threshold.

At the mass limits above there are 388, 157, and 18 parent haloes in the MDLP2/universemachine boxes at redshifts 0.2, 0.75, and 1.25, respectively. These numbers exclude haloes close to the box edges, whose volume as defined in observations (see text below for more details), intersects the box boundaries. The rapid decrease in the number of clusters in each sample is because of the strong evolution of the halo mass function with redshift. These clusters are then projected on to the sky as described in Section 3.3. We choose to place the observer at the same (X, Y) position as the cluster with respect to the reference system of the box. This results in light-cones centred on the each of the selected massive haloes, with a line-of-sight perpendicular to the (X,Y)-plane of the box. We define the FOV in terms of the virial radius of the cluster.

3.4.2 Cluster membership

Although in the simulation box the satellite galaxies associated with a given parent halo are known, we prefer to follow an observational-motivated approach for defining cluster membership based on the projected and radial distances relative to the cluster centre. Cluster member candidates are selected to lie within the projected R_vir radius of the corresponding halo. This is similar to the selection of Martini et al. (2009, 2013). We choose to use R_vir instead of R₂₀₀ since they are similar, but the last is not present in universemachine data set. We checked that this assumption does not affect our final results and conclusions.

The radial distance of clusters member candidates relative to the cluster centre is measured in redshift space as in real observations. Cluster members are those mock galaxies or AGNs with redshift difference to the cluster Δz · c ≤ 3σ_v(1 + z), where z is the redshift of the cluster (fixed to be one of the redshifts of interest, i.e. z = 0.2, 0.75, or 1.25), c represents the speed of light and σ_v is the velocity dispersion of the cluster determined using Munari et al. (2013). This definition corresponds to the selection criteria adopted by Martini et al. (2009) for defining cluster members and it is used for the clusters at redshifts z = 0.2 and 0.75. In the case of Martini et al. (2013) 3σ_v is fixed to |$2000\, {\rm km\,s^{ -1}}$| for all the clusters. This restriction is also adopted in our simulation for the clusters at redshift z = 1.25. This condition defines the volume of the cluster in terms of its velocity dispersion. Fig. 3 shows an example of a simulated observation and demonstrates the impact of selection effects.

Figure 3.

The light-cone of a massive cluster in universemachine at redshift z = 1.25 and halo mass M_vir = 10^14.51M_⊙ h−1. Large circle marks the virial radius of the cluster, all symbols correspond to universemachine galaxies within the light-cone. Different symbols and colours demonstrate selection effects as described in Section 3.4. Purple crosses, red empty squares and blue diamonds indicate galaxies that are excluded from the sample because of the selection effects as indicated in the legend. Green circles and orange triangles indicate the final sample of galaxies and AGNs with |$L_X(2\rm{-}10\, \rm keV)\gt 10^{43}$| erg s⁻¹ after applying all the selection effects.

Open in new tabDownload slide

3.4.3 Galaxy/AGN sample selection

In observations AGNs and/or galaxy samples are typically selected above a given apparent magnitude limit. In Martini et al. (2009, 2013) for example, this is set relative to the knee of the optical and/or mid-infrared luminosity function of galaxies at the corresponding cluster redshift (see Section 2). In simulations, however, like the semi-empirical model described in this work, galaxies are defined by their intrinsic properties, such as stellar mass, SFR, and accretion luminosity. Associating these physical properties to apparent magnitudes requires assumptions on e.g. the SFH of galaxies, the spectral energy distribution (SED) of stellar populations or the shape and normalization of the dust attenuation curve. Assigning SEDs to simulated galaxies is therefore far from trivial and inevitably requires additional modelling steps (e.g. Georgakakis, Ruiz & LaMassa 2020; Pearl et al. 2021)

Our baseline model/observation comparison avoids these additional steps. Instead we make the simplifying assumption that the knee of the observed galaxy optical or mid-infrared luminosity function traces the knee of the underlying galaxy stellar mass function, |$M^{*}_\star$|⁠. This allows us to translate the R band and |$3.6\, \rm \mu$|m apparent magnitude limits of |$0.4\, L^{*}$| adopted by Martini et al. (2009, 2013, see Section 2) to stellar mass cuts. Mock galaxies are selected to be more massive than |$\log \, M^{*}_\star /M_\odot - 0.4\, \rm dex$|⁠. The break of the mass function is fixed to |$\log \, M^{*}_\star =10^{10.7}\, M_\odot$| independent of redshift based on the parametrization of Ilbert et al. (2013). Although the Ilbert et al. (2013) study refers to field galaxies, observations show that the shape of stellar mass function is similar in massive clusters (e.g. Vulcani et al. 2013; Nantais et al. 2016). The translation of the |$0.4\, L^{*}$| apparent magnitude limit to a |$\log \, M^{*}_\star /M_\odot -0.4\, \rm dex$| threshold implies the same average mass-to-light ratio for galaxies. This approximation is justifiable in the case of the high redshift cluster sample (z ≈ 1.25) of Martini et al. (2013), where galaxies are selected in the |$IRAC\, 3.6\, \mu$|m band. At the mean cluster redshift, this wavelength roughly corresponds to rest-fame near-infrared (≈ 1.6 |$\, \mu$|m), where the mass-to-light ratio is not a strong function of the galaxy stellar population. We acknowledge that in the low-redshift (z = 0.2 and 0.75) cluster sample of Martini et al. (2009), galaxies are selected in the R band, where variations of the mass-to-light ratio as a function of the star-formation rate are important. For this sample the approximation of a constant mass-to-light ratio is rough and should be taken with caution.

We further address the limitations above by assigning apparent magnitudes to mock galaxies in the light-cones following the methodology described in Georgakakis et al. (2020). UniverseMachine galaxies are assumed to be described by exponentially declining SFH. The parameters of the SFH model are constrained to reproduce the UniverseMachine stellar masses and instantaneous SFRs of the galaxies at their assigned redshifts in the light-cone. The Bruzual & Charlot (2003) stellar library and the Chabrier (2003) initial mass function are used to synthesize stellar populations for the adopted SFH. The SEDs of star-forming galaxies are extincted by dust assuming the Calzetti et al. (2000) law and E(B − V) = 0.4 mag. The magnitudes of passive galaxies are not extincted by dust. This empirical model is shown to reproduce reasonably well the distribution of apparent magnitudes of galaxies in the COSMOS field (Muzzin et al. 2013).

The assigned apparent magnitudes are sensitive to the instantaneous SFR of mock galaxies. The empirical model of Georgakakis et al. (2020) assumes that star-forming galaxies are on the main sequence of star formation (Schreiber et al. 2015), while quiescent galaxies lie 1 − 1.5 dex below it depending on redshift. This offset for the quenched galaxies is empirically determined to reproduce the observed magnitude distribution of passive galaxies as a function of redshift in the COSMOS field. The universemachine star-forming galaxies are also constrained to follow the main sequence of star formation at different redshifts and are therefore consistent with the assumptions of the empirical SED model of Georgakakis et al. (2020). In contrast, the SFR distribution of quenched galaxies is not as well constrained by observations (see discussion by Pearl et al. 2021). Passive galaxies in universemachine are assigned specific star formation rates (sSFR) that are drawn from a non-evolving lognormal distribution with mean |$\log \rm {sSFR}/\rm {yr}=-11.8$| and scatter 0.36 dex, motivated by observations in the local Universe (Behroozi et al. 2015). This SFR is 2 dex below the main sequence of star-formation at low redshift and therefore inconsistent with the assumptions of the empirical SED model of Georgakakis et al. (2020). This difference in sSFR results in passive galaxies with too faint apparent magnitudes for the empirical model of Georgakakis et al. (2020). In this work, we therefore adopt the definition of quenched galaxies given by equation (2) but then re-scale the corresponding sSFR according to the Georgakakis et al. (2020) prescriptions (i.e. 1–1.5 dex below the main should that be |$M_R^{*}(z=0)$| sequence) before assigning them magnitudes.

Once apparent magnitudes are assigned to mock galaxies in the light-cone, we apply cuts similar to those adopted by Martini et al. (2009, 2013), i.e. one magnitude fainter than the break of the luminosity function in the R (low- and medium-redshift cluster samples, z = 0.2 and 0.75) and 3.6 |$\, \mu$|m (high-redshift cluster samples; z = 1.25) bands. The magnitude limits for the different samples are summarized in Table 1. For the R band in particular, following Martini et al. (2009) we assume that the knee of the luminosity function evolves with redshift as M_R(z) = M_R(z = 0) + z with |$M_R^{*}(z=0)=-21.92$| (Christlein & Zabludoff 2003). This absolute magnitude is converted to apparent magnitude in the R band assuming an elliptical galaxy SED (Ilbert et al. 2009) for the K corrections. The estimated knee of the luminosity function in apparent magnitudes is |$m_R^{*}=17.1$|⁠, 22.3 mag at z = 0.2 and 0.75, respectively. The limits listed in Table 1 are one magnitude fainter than the apparent magnitude of knee of the optical luminosity function at the relevant redshift. For the 3.6 μm band the knee of the luminosity function is assumed to 17.5 mag (apparent magnitude) at z = 1.25 (Mancone et al. 2010). Mock galaxies and AGNs in the simulated light-cones are selected to be brighter than the magnitudes limits listed in Table 1. For the mock AGNs, we further apply the X-ray luminosity limits of Martini et al. (2009, 2013) to define two samples with |$L_X(\rm 2\rm{-}10\, keV) \gt 10^{43}\, \rm {erg \, s^{-1}}$| and |$\rm \gt 10^{44}\, \rm {erg \, s^{-1}}$|⁠.

In the analysis above we ignore the contribution of AGN light to the observed magnitude of galaxies. Modelling this component requires knowledge of the obscuration properties of the AGNs. This latter parameter has the strongest impact on the observed spectral energy distribution of these sources. The specific accretion rate distribution used in this work to seed galaxies with AGN luminosities do not account for the AGN obscuration and therefore cannot be used to model this effect (see also Georgakakis et al. 2020). We simplify this problem by ignoring the AGN contribution to the emission of the galaxy in the optical/mid-infrared. This simplification is equivalent to the assumption that mock AGNs are completely obscured and hence subdominant relative to the stellar emission of galaxies. We will discuss the impact of this assumption on the results and conclusions in the next sections.

4 RESULTS

4.1 Impact of selection effects on the incidence of AGNs

We first explore the sensitivity of the estimate AGN fractions to observational selection effects. Fig. 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 AGNs above a given X-ray luminosity threshold among galaxies more massive than a given stellar mass limit (x-axis of Fig. 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 Fig. 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). Fig. 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 AGNs 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}\, \mathrm{M}_{\odot }$|⁠. This is because in our implementation lower stellar mass hosts require higher λ_sBHAR to produce AGNs 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 Fig. 1, which strongly decreases with increasing λ_sBHAR. For stellar masses |$\gtrsim 10^{11}\, \mathrm{M}_{\odot }$| the duty cycle curves of Fig. 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 Fig. 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.

Figure 4.

Fraction of AGN relative to galaxies above a stellar mass limit as a function of stellar mass. The curves correspond to the predictions of the semi-empirical model described in Section 3.2. Different colours correspond to AGNs more luminous than |$L_X(\rm 2\rm{-}10\, keV)\gt 10^{42}\, erg\, s^{-1}$| (green), |$\rm \gt 10^{43}\, erg\, s^{-1}$| (red) and |$\rm \gt 10^{44}\, erg\, s^{-1}$| (blue). Different line styles indicate different specific accretion rate distribution models adopted to generate the AGN mocks (see Section 3.2). The dashed–dotted lines are for the Georgakakis et al. (2017) and the solid lines correspond to the Aird et al. (2018) specific accretion-rate distributions. Each panels correspond to redshifts from left to right of z = 0.2, 0.75, and 1.25.

Open in new tabDownload slide

For completeness we also show in Fig. 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 AGNs 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 on to the curves of Fig. 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_⋆/M_⊙ > 10.8. This is because this threshold essentially removes a large number of lower mass galaxies, which are typically found in low-mass haloes. 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 Fig. 5.

Figure 5.

Fraction of AGN in galaxies as a function of parent halo mass. The curves correspond to the predictions of the semi-empirical model described in Section 3.2 using the specific accretion rate distribution from Aird et al. (2018). Different colours correspond to AGNs more luminous than |$L_X(\rm 2\rm{-}10\, keV)\gt 10^{42}\, erg\, s^{-1}$| (green), |$\rm \gt 10^{43}\, erg\, s^{-1}$| (red), and |$\rm \gt 10^{44}\, erg\, s^{-1}$| (blue). Different line styles indicate different stellar mass limit cuts M_{⋆, lim} > 10^10.3 M_⊙ (solid line) and M_{⋆, lim} > 10^10.8 M_⊙ (dashed–dotted line). Each panel corresponds to redshifts from left to right of z = 0.2, 0.75, and 1.25.

Open in new tabDownload slide

4.2 X-ray AGN fractions in clusters

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. Fig. 6 plots the fraction of AGN among cluster member galaxies as a function of redshift. The observational results of Martini et al. (2009, 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\rm{-}10\, keV)\gt 10^{43}\, \rm {erg\, s^{-1}}$| and |$\gt 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 realizations 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 per cent confidence interval around the mean. These confidence intervals are the errorbars of the model predictions plotted in Fig. 6. They represent the uncertainty of the mean expected fraction of AGN per cluster. The values of the fractions and its errors for different specific accretion rate models and selection effects are listed in Table 2.

Figure 6.

Evolution of the X-ray AGN fractions in massive clusters. The observations (stars) and semi-empirical model predictions (lines and shaded regions) are plotted at redshifts z = 0.2, 0.75, and 1.25. In both panels the stars (red or blue) are the observationally measured AGN fractions of Martini et al. (2013, see Section 2). Different colours indicate different luminosity cuts. Red is for AGNs with |$L_X (\rm 2\rm{-}10\, keV)\gt 10^{43}\, erg\, s^{-1}$| and blue corresponds to |$L_{X} \rm (2\rm{-}10\, keV)\gt 10^{44}\, erg\, s^{-1}$|⁠. For the shake of clarity the blue stars are shifted by −0.12 in the x-axis direction. The errorbars correspond to the 1σ error. The lower limit at z = 1.25 is due to incompleteness and the upper limit at z = 0.2 corresponds to the 3σ confidence level. In both panels, the lines and corresponding shaded regions represent the predictions of the semi-empirical model on the AGN fraction under different assumptions on the selection effects. On the left-hand panel, the simulated galaxy/AGN samples are selected above the stellar mass limit listed in Table 1. The right-hand panel corresponds to mock AGN/galaxy samples selected above the apparent magnitude thresholds listed in Table 1 (see Section 3.4 for more details). Lines of different colour indicate different luminosity cuts as explained above. The blue lines should therefore be compared with the blue stars and the same for the red lines/symbols. Different line styles indicate different specific accretion rate distribution models adopted to generate the AGN mocks (see Section 3.2). The dashed–dotted lines are for the Georgakakis et al. (2017) and the solid lines correspond to the Aird et al. (2018) specific accretion-rate distributions. The shaded regions within which lines are embedded correspond to the 68 per cent confidence intervals of the mean value calculated using the bootstrapping technique described in the text.

Open in new tabDownload slide

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\rm{-}10\, keV)\gt 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 (∼11 per cent) of the mock clusters have a central galaxy with assigned AGN luminosity |$L_X(\rm 2\rm{-}10\, keV)\gt 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-hand panel in Fig. 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 \gtrsim 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\rm{-}10\, keV)\gt 10^{43}\, \rm {erg\, s^{-1}}$| or |$\gt 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 |$\lesssim 0.1$| dex and therefore does not impact our results and conclusions.

5 DISCUSSION

In this work, we study the fraction of AGN in massive clusters with a semi-empirical modelling technique that allows the generation of realistic AGN mock catalogues. Dark matter haloes from cosmological simulations are seeded with galaxies using abundance matching techniques (Behroozi et al. 2019). On top of these galaxies accretion events by supermassive black holes are painted. The latter step is based on state-of-the-art specific accretion rate distributions derived from observations (Georgakakis et al. 2017; Aird et al. 2018). The zero-order assumption of the semi-empirical model is that the incidence of accretion events in galaxies is independent of environment, i.e. galaxies are seeded with AGNs using the same empirical relations independent of the mass of the parent halo. We refer to the predictions of the model as ‘field expectation’ for the reasons explained in Section 3.2. This methodology reproduces by construction the halo mass function, stellar mass function and the X-ray luminosity function as demonstrated in Fig. 2.

5.1 AGN fractions in high-redshift clusters

A striking result in Fig. 6 is the higher observed fraction of AGNs in massive clusters at z = 1.25 compared to the model predictions. The largest discrepancy of nearly 1 dex is for powerful AGNs with |$L_X(\rm 2\rm{-}10\, keV) \gt 10^{44} \, erg \, s^{-1}$|⁠. An enhanced fraction is also observed for moderate luminosity AGN, |$L_X(\rm 2\rm{-}10\, keV) \gt 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 ≈ 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 Fig. 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 Fig. 7 are normalized 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 haloes) and by construction reproduces the observed stellar mass function at different redshifts estimated using extragalactic survey fields (e.g. Ilbert et al. 2013; Moustakas et al. 2013; Muzzin 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 Fig. 6 cannot be attributed to systematic variations of the total stellar mass function with environment.

Figure 7.

Stellar mass function normalized to the total stellar mass of the sample. Each panel correspond to one of the redshifts z = 0.2, 0.75, and 1.25. The blue curve is for all galaxies in the corresponding universemachine box and therefore represents the field stellar mass function. The orange colour is the stellar mass function of the satellites of the selected clusters (see Section 3.4 for more details). The shaded orange region correspond to the Poisson noise uncertainty.

Open in new tabDownload slide

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 Fig. 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\gtrsim 1$| clusters compared to the field expectation in Fig. 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 centred 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 ≈ 1. They report a statistically significant overdensity of X-ray selected AGNs within the infall cluster region at a project distance of |$\sim 2\rm{-}2.5 \, R_{500}$| relative to the cluster centre. 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\gtrsim 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) identify AGN by fitting model templates to the multiwavelength SEDs of mid-infrared selected cluster members.

5.2 Reproducing high-z AGN fractions in clusters

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 ≈ 1.25 shown in Fig. 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 ≈ 1.25 shown in Fig. 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\gt \rm 10^{44}\, erg \, s^{-1}$| data point or larger than the lower limit for |$L_X\gt \rm 10^{43}\, erg \, s^{-1}$| AGN in Fig. 6. The general trend is that broader distributions (larger scatter) require lower means λ_sBHAR to reproduce the data. We also find that Gaussians with mean of |$\lambda _{sBHAR} \gtrsim 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}\lesssim 10^{-3}$| produce too few luminous AGN for any scatter value and are also rejected.

Fig. 8 (left-hand 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}\lesssim \lambda _{sBHAR} \lesssim 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 Fig. 8 make different predictions on the number of cluster AGN as a function of accretion luminosity. This is demonstrated in the middle panel of Fig. 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 normalization from the field one, also plotted in Fig. 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 normalization XLFs shown in Fig. 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.

Figure 8.

Comparison of Gaussian and the observationally derived (see Section 3.2) specific accretion rate distribution models for the cluster AGN sample described in Section 3.4 at z = 1.25. The left-hand panel shows the specific accretion rate distributions that describe the probability of a galaxy hosting an AGN with specific accretion rate λ_sBHAR at 1 < z < 1.5 and stellar masses 10^10.5 < M_⋆/M_⊙ < 10¹¹. The observationally derived models are the purple (Aird et al. 2018) and green (Georgakakis et al. 2017) curves. Shaded region indicates 68 per cent confidence intervals. The yellow, light brown, and dark brown curves correspond to the Gaussian distributions (see text for details) with different means and dispersions as indicated in the legend. The middle panel compares the X-ray luminosity function for the cluster sample, predicted by the specific accretion rate distribution models (lines) with observations (orange points; Georgakakis et al. 2017). The different lines indicate the observationally derived (representative of the field population) specific accretion rate models (solid, Georgakakis et al. 2017; Aird et al. 2018, green and purple, respectively,) and the Gaussian specific accretion-rate models, proposed to match the AGN fraction in clusters at this redshift (dashed, with different colors indicating different Gaussian parameters as indicated in the legend). The right-hand panel shows the predicted AGN fraction in clusters at z = 1.25 following the selection effects as explained in Section 3.4. The stars (red or blue) are the observationally measured AGN fractions of Martini et al. (2009, 2013). Different colours indicate different luminosity cut. Red is for AGN with |$L_X (\rm 2\rm{-}10\, keV)\gt 10^{43}\, erg\, s^{-1}$| and blue corresponds to |$L_{X} \rm (2\rm{-}10\, keV)\gt 10^{44}\, erg\, s^{-1}$|⁠. The predictions from the observationally derived specicific accretion-rate models in Fig. 6 are shown here with horizontal dashed–dotted (Georgakakis et al. 2017) and solid lines (Aird et al. 2018). The circle, diamond, and plus markers correspond to the predictions of the three Gaussian specific accretion-rate models shown in the left-hand panel with parameters as indicated in the legend.

Open in new tabDownload slide

Finally, the difference between field and cluster specific accretion rates plotted in the left-hand panel of Fig. 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 on to the supermassive black hole.

5.3 Evolution of AGN incidence in cluster versus field

Contrary to the results at z = 1.25 discussed above, Fig. 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 AGNs or the type of observational selection effects (mass versus 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 \rm{-} 10 \, keV) \gt 10^{43} \, erg \, s^{-1}$| predicted by the model compared to the observations. For more luminous AGN, |$L_X (\rm 2 \rm{-} 10 \, keV) \gt 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 flavour a large fraction (≈10 per cent; 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 (≈1–3 per cent) 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 \rm{-} 10 \, keV) \gt 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\gtrsim 1$| are found to be more efficient in triggering accretion events on to supermassive black holes compared to less dense regions. Such an environmental dependence is not present at z ≈ 0.75 and possibly inverses at low redshift, z ≈ 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 ≈ 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 protoclusters at higher redshift |$z\gtrsim 2$| find that AGNs 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 \lt 1\, Mpc$|⁠, in agreement with our findings.

It is also interesting that observational studies of the Halo Occupation Distribution (HOD) of AGNs are broadly consistent with this picture. At low-redshift, |$z\lesssim 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.

5.4 Physical interpretation

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 normalized to halo mass increases with look-back 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 versus density relation (Butcher & Oemler 1984), whereby cluster cores at |$z\gtrsim 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\gtrsim 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\gtrsim 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. Vulcani et al. 2018; Kotecha et al. 2022) 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₅₀₀) of massive clusters (⁠|$M_{500} \gt 10^{14}\, M_\odot$|⁠) at z ≈ 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 haloes. 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:

• 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 haloes as well as the accretion rate on to them increases with redshift (e.g. Gottlöber, Klypin & Kravtsov 2001; Fakhouri & Ma 2008; McBride, Fakhouri & Ma 2009). Therefore the flow of galaxies from the cosmic web on to massive haloes is expected to be higher at earlier times therefore leading to enhanced galaxy interaction rates. The semi-analytical model of Gatti et al. (2016) suggests that galaxy interactions dominate relatively luminous AGN |$L_X\gtrsim 10^{43} \rm \, erg \, s^{-1}$| in massive haloes (⁠|$\gtrsim 10^{14}\, M_\odot$|⁠) and low/intermediate redshift, |$z\lesssim 1$|⁠. Instead, internal processes, such as disc 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 disc instabilities dominating in massive haloes. This is because of the higher gas and disc fraction of galaxies at these higher redshift (Gatti et al. 2016).

• Ram pressure. Interaction between the intracluster 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 on to 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\rm{-}51{{\ \rm per\ cent}}$|⁠) 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 favourable 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 ≈ 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 haloes 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\gtrsim 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 on to massive haloes 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 haloes are typically less massive. Such an evolution pattern is consistent with the higher fraction of AGN in massive clusters at |$z\gtrsim 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\gtrsim 1$|⁠.

6 SUMMARY AND CONCLUSIONS

This paper addresses the fundamental question of the role of small-scale environment (⁠|$\rm \lt 1\, Mpc$|⁠) in triggering accretion events on to the supermassive black holes at the nuclear regions of galaxies. We tackle this issue by developing a flexible semi-empirical model that populates the dark matter haloes of cosmological N-body simulations with AGN using observational relations on the incidence of accretion events among galaxies. This zero-order assumption of the model is that the probability of accretion events are independent of the parent halo masses of galaxies, i.e. agnostic to their small-scale environment. Moreover, the observationally derived AGN incidence probabilities adopted by the model are representative of the field galaxy population, i.e. those outside massive groups or clusters. This model is used to predict the fraction of AGN in massive clusters of galaxies at different redshifts and compare against observational results by carefully taking into account observational selections effects. Any differences between model predictions and observations would point to an environmental dependence of AGN triggering mechanisms. Our main findings are

• the X-ray AGN fraction in massive clusters are larger than the model prediction at z ∼ 1.25. This points to a strong environmental dependence of AGN triggering at high redshift. Black hole accretion events are promoted in massive haloes relative to the model expectation, which in turn represents the field expectation.

• the model predictions are consistent with the observed AGN fractions of interemediate redshift (z ∼ 0.75) clusters of galaxies. At this redshift, it appears that massive haloes do not promote or suppress AGN activity relative to the field predictions of the model.

• at low redshift, z ∼ 0.2, the model overpredicts the fraction of |$L_X(\rm 2\rm{-}10\, keV) \gt 10^{43} \, erg \, s^{-1}$| AGN in clusters compared to observations. This suggests a suppression of AGN activity in clusters relative to the field expectation of the model at z ∼ 0.2.

• overall the points above suggest a differential redshift evolution of the AGN fraction in clusters relative to the field predictions of our semi-empirical model.

The observed trends above may be related to the increasing gas content of galaxies with increasing redshift, coupled with mechanisms such as galaxy interactions or ram-pressure. Both of these processes under certain conditions could promote AGN activity among galaxies that fall on to massive clusters at higher redshift.

ACKNOWLEDGEMENTS

The authors would like to thank the anonymous referee for their comments to the paper, Ángel Ruiz Camuñas for discussions that significantly improved the simulation code used in this work and Peter Behroozi for his help with the UniverseMachine dataset. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement no 860744. AL is partly supported by the PRIN MIUR 2017 prot. 20173ML3WW 002 ‘Opening the ALMA window on the cosmic evolution of gas, stars, and massive black holes’. SB acknowledges the project PGC2018-097585-B-C22, MINECO/FEDER, UE of the Spanish Ministerio de Economia, Industria y Competitividad. This research made use of Astropy,⁴ a community-developed core python package for Astronomy (Astropy Collaboration 2013, 2018); colossus⁵ (Diemer 2018), numpy⁶ (van der Walt, Colbert & Varoquaux 2011), scipy⁷ (Virtanen et al. 2020), Matplotlib⁸ (Hunter 2007) and HaloMod⁹ (Murray, Power & Robotham 2013; Murray et al. 2021). For the purpose of open access, the authors have applied a CC-BY public copyright license to any author accepted manuscript version arising.

DATA AVAILABILITY

The data products and relevant code to reproduce the results of this paper are available at https://github.com/IvanMuro/agn_frac_data_release.

Footnotes

REFERENCES