Scrutinizing evidence for the triggering of active galactic nuclei in the outskirts of massive galaxy clusters at z ≈ 1

ABSTRACT

Environmental effects are believed to play an important yet poorly understood role in triggering accretion events onto the supermassive black holes (SMBHs) of galaxies (active galactic nuclei; AGNs). Massive clusters, which represent the densest structures in the Universe, provide an excellent laboratory to isolate environmental effects and study their impact on black hole growth. In this work, we critically review observational evidence for the preferential activation of SMBHs in the outskirts of galaxy clusters. We develop a semi-empirical model under the assumption that the incidence of AGN in galaxies is independent of environment. We demonstrate that the model is broadly consistent with recent observations on the AGN halo occupation at z =

0.2, although it may overpredict satellite AGN in massive haloes at that low redshift. We then use this model to interpret the projected radial distribution of X-ray sources around high redshift (z ≈ 1) massive (⁠|$\gt 5 \times 10^{14} \, M_\odot$|⁠) clusters, which show excess counts outside their virial radius. Such an excess naturally arises in our model as a result of sample variance. Up to 20  per cent of the simulated projected radial distributions show excess counts similar to the observations, which are however, because of background/foreground AGN and hence, not physically associated with the cluster. Our analysis emphasizes the importance of projection effects and shows that current observations of z ≈ 1 clusters remain inconclusive on the activation of SMBHs during infall.

1 INTRODUCTION

A major challenge in current astrophysical research is to understand the formation and evolution of galaxies in the Universe.The difficulty in addressing this issue is that the relevant physical processes, such as the cooling of gas, the formation of stars, and the injection of energy and metals into the interstellar medium by e.g. dying stars, are complex, interconnected, and operate over a wide range of temporal and spatial scales (e.g. Benson 2010; Somerville & Davé 2015). A development that has changed the way we view galaxy evolution has been the realization that nearly every spheroidal galaxy hosts at its nuclear regions a black hole with a mass that may exceed |$10^9 \, \rm {M}_\odot$| that appears to correlate with the mass of the stellar population (e.g. Graham et al. 2011; Kormendy & Ho 2013). Although these supermassive black holes (SMBHs) are gravitationally insignificant for their galaxies, theoretical arguments and observational results suggest that their energy output during their growth phases has a strong impact on the interstellar medium and can affect the evolution path of their hosts (e.g. Fabian 2012). Understanding in detail this symbiotic relationship is therefore important for painting a complete picture of galaxy evolution. A first step toward this goal is to understand the physical conditions that produce accretion flows onto SMBHs, thereby leading to their growth and the release of energy that is observed as active galactic nuclei (AGNs).

The activation of SMBHs relies on two factors. The availability of cold gas in the galaxy to fuel these compact objects and a mechanism that is able to drive this material to the galactic centre in the vicinity of the SMBH. Secular processes that occur during the lifetime of galaxies can generate conditions that fulfil the requirements above and hence promote the growth of black holes. For example, recycled gas produced by normal stellar evolution can provide sufficient reservoirs of available fuel and lead to recursive cycles of black hole accretion flows (e.g. Ciotti & Ostriker 2007). Disc instabilities (e.g. Hopkins & Hernquist 2006; Gatti et al. 2016) and galactic bars (Cisternas et al. 2015) are efficient in removing angular momentum from the interstellar gas thereby driving it towards the central regions of the galaxy where it can be accreted by the SMBH. In the early Universe, the direct collapse of low angular momentum gaseous baryons is proposed to lead to starburst events as well as the rapid growth of black holes and ultimately produce the progenitors of present-day early-type massive galaxies (Shi et al. 2017; Lapi et al. 2018). In addition to the in-situ processes above, environmental effects are also thought to play an important role in modulating accretion flows onto SMBHs. For example, in the case of galaxies in dense regions of the cosmic web, ram pressure may initially compress the cold gas in the nuclear regions of galaxies (Marshall et al. 2017; Ricarte et al. 2020) and hence, promote accretion flows onto SMBHs (Peluso et al. 2022). In the longer term, however, this process acts to deplete the cold gas reservoirs of the galaxies (Steinhauser, Schindler & Springel 2016), thereby suppressing the growth of their black holes. Gravitational interactions and mergers are long thought to represent an important AGN trigger (e.g. Di Matteo, Springel & Hernquist 2005; Koulouridis et al. 2006, 2013; Gatti et al. 2016) and perhaps the dominant mechanism in the case of the most luminous SMBH accretion events (e.g. Glikman et al. 2015; Araujo et al. 2023). At high redshift flows of cold gas from the cosmic web onto galaxies are proposed to be common leading to both intense star-formation and AGN activity (Bournaud et al. 2012; DeGraf et al. 2016).

This paper focuses on AGN triggering mechanisms that are pertinent to the densest structures in the Universe, massive galaxy clusters. These systems offer a perfect laboratory for isolating environmental effects to explore how they modulate black hole growth. Additionally, by scanning clusters of galaxies from beyond their outskirts to their cores, it is possible to sample a broad range of densities and therefore witness the onset of environmental effects and study their impact as a function of local density. Most observational evidence of the fraction of AGN in low redshift (⁠|$z\lesssim 0.3$|⁠) clusters indicate that the nuclear activity in galaxies is suppressed in these dense environments, particularly close to the centre of the potential well (Haines et al. 2012; Martini et al. 2013; Sabater, Best & Argudo-Fernández 2013; Koulouridis et al. 2014; Mishra & Dai 2020). Some observations find no or little impact of the environment on AGN activity at low redshift (e.g. Miller et al. 2003; Haggard et al. 2010; Pimbblet et al. 2013). However, at least part of the discrepancy can be attributed to the different selection effects (e.g. X-ray selection versus optical). At intermediate redshifts (z ≈ 0.7), the difference between cluster environment and field seems to be level off (see e.g. Eastman et al. 2007; Galametz et al. 2009; Martini, Sivakoff & Mulchaey 2009; Martini et al. 2013; Ehlert et al. 2014; Bufanda et al. 2017) and perhaps reverse at |$z\gtrsim$|1 (see Lehmer et al. 2009; Digby-North et al. 2010; Krishnan et al. 2017; Tozzi et al. 2022; Monson et al. 2023; Muñoz Rodríguez et al. 2023; Toba et al. 2024). These studies show that the incidence of AGN in clusters of galaxies is similar to the field expectation at z ≈ 0.3–0.8 and exceeds this value at earlier cosmic times. The cluster environment, therefore, appears to promote black hole growth outside the local Universe. Efficient activation of AGN in dense regions points to physical mechanisms that operate preferentially in these environments, such as a higher galaxy interaction rate (Gatti et al. 2016) or ram-pressure (Poggianti et al. 2017; Peluso et al. 2022). These processes are expected to be more efficient at the outskirts of clusters (e.g. Toba et al. 2024) where the local density is lower and the relative velocities of galaxies smaller.

An infalling population of active black holes may imprint observable features on the radial distribution of AGN within a cluster (e.g. Rihtaršič et al. 2024). There is indeed evidence that the fraction of AGN relative to galaxies is decreasing towards the cluster centre suggesting a higher incidence of AGN at the cluster outskirts (Martini et al. 2009; Ehlert et al. 2014; de Souza et al. 2016; Lopes, Ribeiro & Rembold 2017; Mishra & Dai 2020; Stroe & Sobral 2021; Koulouridis, Gkini & Drigga 2024). Additionally, there are claims that the projected counts of AGN show an excess outside the virial radius of clusters (Johnson, Best & Almaini 2003; Ruderman & Ebeling 2005; Fassbender, Šuhada & Nastasi 2012; Koulouridis & Bartalucci 2019), which could be interpreted as direct evidence of SMBH activation during infall. However, these results remain controversial with a number of studies failing to observe such projected overdensities (Ehlert et al. 2014; Mo et al. 2018; Mishra & Dai 2020). Part of the discrepancy can be attributed to differences in cluster halo masses or cluster dynamical states among the various samples (e.g. Hashiguchi et al. 2023), AGN selection effects such as flux limits or selection wavelength as well as cluster to cluster variations (see Martini, Mulchaey & Kelson 2007).

In this work, we revisit claims for an excess of AGN activity in the outskirts of clusters by developing a semi-empirical model to interpret the observed X-ray AGN radial distributions presented by Koulouridis & Bartalucci (2019). This work uses Chandra observations of a well-defined sample of clusters with carefully measured masses and sizes to find a statistically significant excess of X-ray point sources at a distance of about |$2.5\, R_{500}$| from the cluster centre, where R₅₀₀ is the radius that encloses a volume with mass density 500 times the critical one of the Universe at the redshift of interest. The Koulouridis & Bartalucci (2019) work has a number of key features that greatly facilitate the modelling and interpretation. The first is the transparent selection of the clusters and the corresponding AGN which can be replicated in the modelling. The second is the fact that the radial distributions are expressed in units of R₅₀₀, thereby allowing direct comparison of clusters with different masses and extents. The semi-empirical modelling approach we develop in this paper provides an excellent handle on systematics and selection effects and enables us to explore the impact of projection effects and sample variance in the radial distributions of AGN in clusters of galaxies. Our modelling is based on observationally derived relations to populate dark matter haloes extracted from N-body simulations with AGN and galaxies under the assumption that the incidence of AGN does not depend on environment. The comparison with the observations follows the principles of forward modelling to generate realistic cluster observations that include selection effects such as flux limits and the finite Chandra field of view. Section 2 presents the observations and the cluster sample used in this work. 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 the current debate on AGN radial distribution in clusters. We adopt a flat ΛCDM cosmology with parameters Ω_m = 0.307, |$\Omega _\Lambda$|= 0.693, h = 0.678 consistent with the Planck results (Planck Collaboration XIII 2016).

2 OBSERVATIONS

This work uses Chandra X-ray observations of massive clusters of galaxies at z ≈ 1 presented by Koulouridis & Bartalucci (2019). Their sample is selected using the Sunyaev–Zeldovich effect (SZ) and it is composed by the five most massive clusters (⁠|$M_{500}^{SZ}\gtrsim 5\times 10^{14}\rm {M}_\odot$|⁠) in the South Pole Telescope and Planck catalogues at that redshift (see Bleem et al. 2015; Planck Collaboration XIII 2016). These are the only clusters at that redshift for which detailed analysis of their X-ray profiles have been carried out (Bartalucci et al. 2017) to provide robust constraints on their R₅₀₀ (0.7–1 Mpc) and M₅₀₀ (mass range |$3-8 \times 10^{14}\, M_\odot$|⁠). Koulouridis & Bartalucci (2019) explore the projected radial distribution of X-ray sources in their cluster sample and find evidence for a systematic excess of counts at a projected radius of |$\approx 2.5\, R_{500}$|⁠. In this paper, we build a semi-empirical model of the radial distribution of AGN in massive dark matter haloes identified in cosmological N-body simulations. Then we compare the predictions of the model with the observational results of Koulouridis & Bartalucci (2019), using the principles of forward modelling. In this approach, observational effects such as X-ray flux limits and the Chandra field of view are included in the modelling to generate simulated data sets that mimic real observations. In that respect, an important part of the simulations is the X-ray selection function of the observations, which measures the probability of detecting X-ray point sources of a given flux as a function of position within the Chandra footprint. For that reason we choose to re-analyse the Chandra X-ray observations used by Koulouridis & Bartalucci (2019) with the reduction pipeline presented by Laird et al. (2009) and Nandra et al. (2015). The key feature of this pipeline is the sensitivity maps that are constructed following methods presented in Georgakakis et al. (2008) and quantify to a high level of accuracy the selection function of the detected X-ray sources.

In brief, the reduction uses standard CIAO tasks to analyse the raw Chandra/ACIS-I imaging data and produce level-2 event files for individual pointings. Multiple observations of the same field are merged to generate a single event file as well as co-added images and exposure maps in four energy bands 0.5–7.0 keV (full), 0.5–2.0 keV (soft), 2.0–7.0 keV (hard), and 4.0–7.0 keV (ultrahard). Sources are detected independently in each of these spectral intervals following a two-pass process. A seed catalogue of candidate sources is first constructed using the CIAO wavelet-based source detection task WAVDETECT at a low detection threshold of 10⁻⁴. This is to ensure a high level of completeness of the source list. Photons (source and background) at the position of each candidate source are then extracted within apertures that correspond to the 70 per cent encircled energy fraction (EEF) radius of the Chandra point spread function (PSF) at the source position. The expected background level in each aperture and spectral band is also measured after removing the contribution of nearby source photons. Finally, we estimate for each source the Poisson probability that the observed number of photons within the aperture is the result of background fluctuations. The final catalogue in a given spectral band contains those X-ray sources with Poisson probability as defined above <4 × 10⁻⁶. The choices of the EEF radius for aperture photometry and the Poisson probability cutoff are a trade-off between completeness and purity (e.g. Laird et al. 2009). Larger values of EEF radii would reduce the completeness of the catalogue because faint sources may be swamped within the higher level of background of the bigger apertures. Smaller values of EEF, on the other hand, would increase the number of spurious sources (at fixed Poisson probability threshold). Reducing the Poisson probability cutoff renders the source detection algorithm more conservative with higher purity but also lower completeness and vice versa.

X-ray fluxes are determined assuming a power-law X-ray spectrum with photon index Γ = 1.4 (similar to the diffuse X-ray background; Akylas et al. 2012) absorbed by the Galactic hydrogen column density appropriate for each field. The pipeline also produces sensitivity maps (see Georgakakis et al. 2008), which measure the probability of detecting an X-ray source with a given count rate or flux as a function of position within the surveyed area. In this work, we use the 1D representation of the sensitivity map, the X-ray area curve, which provides an estimate of the total survey area in which a source with a given count rate or flux can be detected.

Next, we describe the construction of the radial distribution of X-ray sources in each of the clusters in the sample of Koulouridis & Bartalucci (2019). We use X-ray sources selected in the full-band (0.5–7 keV) and group them in radial bins of width 0.5  · R₅₀₀. These radii are estimated by Bartalucci et al. (2018) by modelling the extended X-ray emission profile of each cluster. The determination of the projected radial distribution of X-ray sources requires the statistical subtraction of the expected number density of foreground/background X-ray sources. This is determined using the number counts as a function of the flux of the extragalactic field X-ray source population, i.e. their log N − log S distribution. For a given cluster and R₅₀₀ radial bin, i, we first determine the full-band sensitivity curve of the ring with inner and outer radius of i/2 · R₅₀₀ and (i + 1)/2 · R₅₀₀ following the methods described in Georgakakis et al. (2008). We then convolve this with the differential full-band number counts presented by Georgakakis et al. (2008). This calculation yields the number of extragalactic field X-ray sources (i.e. not associated with the cluster) expected to be detected in the ring under consideration at the depth of the specific Chandra observation. This expectation value is then subtracted from the observed number of X-ray sources in the ring. The resulting distributions for the clusters PLCKG266.6-27.3 and SPTCLJ2146-4633 are shown in Fig. 1. The two selected clusters are the ones that show the highest overdensity at |$2.5\, R_{500}$| in the sample of Koulouridis & Bartalucci (2019). This figure also shows that our re-analysis confirms the results of Koulouridis & Bartalucci (2019).

In the next sections, we will use the full-band sensitivity maps of the clusters PLCKG266.6-27.3 (⁠|$\rm {M}_{500}=8.38^{+0.35}_{-0.36}~\rm {M}_\odot$|⁠, |${R}_{500}=993\pm 14\, \rm {kpc}$|⁠, and z = 0.972; Bartalucci et al. 2018) and SPTCLJ2146-4633 (⁠|$\rm {M}_{500}=3.15^{+0.13}_{-0.14}~\rm {M}_\odot$|⁠, |${R}_{500}=728^{+10}_{-11}\, \rm {kpc}$|⁠; and z = 0.933; Bartalucci et al. 2018) to forward model the X-ray selection function of real Chandra observations. The first represents a deep, ≈200 ks, X-ray observation, the second corresponds to a shallower Chandra data set, ≈70 ks. We will use the sensitivity maps of these observations to explore the impact of different X-ray depths on our results and conclusions. We reiterate that we choose these two clusters because they are the one in the sample of Koulouridis & Bartalucci (2019) that show the largest amplitude excess counts in their projected radial distributions, which are interpreted as evidence for AGN triggering in their outskirts.

3 METHODOLOGY

3.1 The semi-empirical model of AGN and galaxies

In this section, we describe the development of the semi-empirical model that is used to interpret the observations presented in Section 2 on the radial distribution of AGN in massive clusters of galaxies. The semi-empirical approach is a flexible data-driven method that produces realistic mock catalogues of galaxies (e.g. Moster, Naab & White 2018; Behroozi et al. 2019; Grylls et al. 2019) and/or AGN (e.g. Comparat et al. 2019, 2020; Seppi et al. 2022; Zhang et al. 2022). By construction, such mocks obey observed properties of galaxy and/or AGN populations, e.g. the stellar mass function, the star formation main sequence at different cosmic times, and the AGN luminosity function. In contrast with other modelling methods, such as hydrodynamical simulations or semi-analytical models, the semi-empirical approach does not rely on a set of recipes to describe the physical mechanisms that regulate galaxy/AGN evolution. Instead, empirical assumptions are made, e.g. the stellar mass of a galaxy correlates with halo mass, which can usually be described by few parameters. Because of its simplicity, the semi-empirical approach is ideal for testing specific hypotheses by comparing simulations with observations. It is this latter point that motivates the use of the semi-empirical approach in our analysis, instead of more complex and physically driven simulations of massive clusters of galaxies (e.g. Cui et al. 2018).

In this work, we follow the methodology described in Muñoz Rodríguez et al. (2023) to construct AGN mock catalogues. In brief, the backbone of the model is a dark matter only N-body simulation. It provides the dark matter halo structure within which galaxies form and evolve. We choose to use the MultiDark PLanck2 (MDPL2; Klypin et al. 2016) because it is one of the largest volume, high resolution, and public cosmological simulations. It has a box size of 1000 Mpc h⁻¹, a mass resolution of |$1.5\times 10^9\,$|M_⊙ h⁻¹, and 3 840³ (∼57 × 10⁹) particles. Dark matter haloes are populated with galaxies using abundance matching techniques. In particular, we use the universemachine model of Behroozi et al. (2019) implemented for the MDPL2 dark matter N-body simulation. This model assigns galaxies to haloes by parametrizing the star-formation rate (SFR) in terms of halo properties (mass, accretion history, and cosmic time). By integrating the SFR across the halo history, it is possible to predict observables that are compared with real observations, including for example, the stellar mass function and the evolution of the cosmic star-formation rate density. The best model is found by iterating the comparison between predictions and observations to explore the model parameter space. The end product of universemachine are catalogues of dark matter haloes, each of which is assigned a galaxy stellar mass and a SFR.

Following Muñoz Rodríguez et al. (2023), an AGN luminosity is assigned to each galaxy in universemachine using observational measurements of the AGN specific accretion rate distribution (SARD; Aird et al. 2012; Bongiorno et al. 2012, 2016; Georgakakis et al. 2017; Aird, Coil & Georgakakis 2018). This quantity describes the probability of a galaxy hosting an accretion event onto its supermassive black hole with specific accretion rate (SAR) λ_SAR = L_X/M_⋆, where L_X is the AGN luminosity (in this case at X-rays) and M_⋆ is the stellar mass of the host galaxy. The observationally derived SARDs are used to assign accretion events to mock galaxies in a probabilistic way (e.g. Georgakakis et al. 2019; Aird & Coil 2021; Muñoz Rodríguez et al. 2023) and therefore include mock AGN in the universemachine boxes. The fundamental assumption of this step is the lack of a physical connection between the accretion events and the environment, i.e. the AGN incidence is stochastic in nature and independent of the halo mass. The process of constructing the galaxy and AGN semi-empirical model described above is illustrated in the first three panels from left to right on the upper branch of Fig. 2.

The catalogues of mock AGN and galaxies produced above need to be further processed to mimic observations of the real Universe and allow the comparison with observational results in a forward modelling manner. The essential step for achieving this is the projection of the boxes onto the sky plane to construct light cones as in Muñoz Rodríguez et al. (2023). However, the light-cone requirements of this work are very different from those in Muñoz Rodríguez et al. (2023). As a result the construction of this product deviates from our previous study and is described in detail in the next section.

3.2 Light cones

In this work, we explore the projected radial distribution of X-ray-selected AGN in galaxy clusters and how this is affected by sample variance. At the simulation level, this is investigated by generating light cones of massive dark matter haloes whose sightlines probe different paths through the cosmic web. This is demonstrated in Fig. 3 which shows two different sightlines to a particular halo (left panel). The corresponding projected structures along these sightlines are also shown in the figure. In the next sections, we first discuss the general approach for constructing light cones (Section 3.2.1) and explain how this is modified to allow more freedom in the choice of sightlines to a particular halo (Section 3.2.2).

3.2.1 General light-cone construction

Extragalactic surveys are typically characterized by a finite field of view and a flux limit at some waveband that allows the detection of astrophysical sources (galaxies or AGN) over a wide range of redshifts. Dark matter N-body simulations like MDPL2 have a finite box size, which when projected onto the sky plane samples only a limited redshift range.¹ Producing mocks over a wide redshift interval requires that simulated boxes are used as building blocks to construct a 3D pavement. The stack of boxes can be extended from z = 0 to an arbitrary maximum redshift (z_max) by selecting an appropriate number of boxes. A wide range of redshift, however, corresponds to a significant look-back time, during which the structure of the Universe evolves strongly. This effect can be captured by selecting different dark matter simulation boxes that correspond to different redshifts. They represent snapshots of the cosmic web at distinct times during the lifetime of the Universe. Using different snapshots and stacking them we construct catalogues that describe the evolution of the structure in the Universe over a wide range of redshifts. The skeleton of these catalogues can be described as an onion-shell structure where each slice corresponds to a different snapshot. A potential issue with this approach is that since distinct snapshots represent the same volume of the simulated universe, the positions of specific structures are correlated between different boxes. This is known as the repetition problem. Diverse alternatives have been proposed in the literature to address this limitation. We implement random tiling, which decorrelates relative positions between different boxes by rotating them along the main axes of the box when they are stacked (see Blaizot et al. 2005; Bernyk et al. 2016). This process introduces spurious correlations at scales bigger than each redshift slice, however, note that this is a second-order effect for our particular application and has little or no impact on the sample variance.

The origin of the reference system of a given simulation box is assumed to be located at the centre of the box. The Cartesian coordinates of the individual objects therefore, take values in the range [−L_box/2, L_box/2], where L_box is the length side of the simulated volume. The hypothetical observer is located onto the XY plane at Z = 0. Its precise position on the plane can be almost arbitrary. Constraints are discussed in Section 3.3. The box is located at a comoving distance that corresponds to a reference redshift (z_ref) with respect to the position of the observer. This is achieved by offsetting the box along the Z-axis by Δz defined as

where D_c(z = z_ref) is the comoving distance that corresponds to z_ref. The coordinates in the equation above are defined as x = x_pp − x_RS, y = y_pp − y_RS, and z = z_pp, where RS are the coordinates of the observer and PP are the coordinates of the pivot point. The latter are defined as the coordinates of a point within the box that has a distance of exactly D_c(z = z_ref) from the observer. The exact location of the pivot point within the box is a free parameter although it is typically chosen to be the centre of the box. The z_ref usually corresponds to the reference redshift of the snapshot. Deviations from these norms are discussed in Section 3.3.

The stacking of boxes requires some overlap between consecutive boxes to avoid empty volumes which would generate an incomplete light cone. This is achieved by imposing the condition D_c(z_i) − D_c(z_{i + 1}) < L_box, i.e. the comoving distance between the centres of consecutive boxes should be smaller than their comoving length. However, the overlap produces artificial overdense regions because the same volume contains objects from two different boxes. This is demonstrated on the left panel of Fig. 4, which shows the stack of two boxes. The intersecting volume contains the individual haloes of each box and therefore, it has an artificially enhanced density. We address this issue by defining a boundary surface of constant comoving distance (or redshift) relative to the observer. We refer to this as the split-overlap-surface (sos). It determines which objects are adopted from each box. Above the surface, only haloes from the box on the top are kept. Whereas below the surface, only objects from the bottom box are retained. This is illustrated in the middle panel of Fig. 4, where the lower curved line represents the split-overlap-surface.

The split-overlap-surfaces define a set of box slices, i.e. the onion shell structure of the light cone. The stacking of the simulation boxes to construct the light cone follows a top-to-bottom approach: the pivot point of the highest redshift box is defined and the appropriate offset relative to the observer is applied to it. The sightline between that pivot point and the observer define the axis of symmetry of the light cone. Lower redshift boxes are then added underneath the first one by defining appropriate pivot points and split-overlap-surfaces. The relative angle between the objects in the box slices and the selected sightline is calculated. This angle can be decomposed into a right ascension and declination on the unit sphere. The redshifts associated with the individual haloes correspond to their comoving distances with respect to the observer. Then the input field of view of the light cone is applied by rejecting objects with angular distances larger than the adopted solid angle. This is illustrated on the right panel of Fig. 4, where only objects within the limits of the field of view are included.

3.2.2 Cut-and-paste method

For our specific application, it is necessary to construct light cones that intersect a particular halo position (i.e. that of a massive cluster) at a comoving distance from the fiducial observer that corresponds to a fixed redshift. Therefore the pivot point of the box that contains this particular halo is set to the Cartesian position of this halo. In this case, the methodology described in the previous section has a limitation that is demonstrated in the left panel of Fig. 5. The sightline to the target object may intersect the boundaries of a box in the stack before reaching the maximum redshift of the light cone. We address this issue by modifying the methodology described in Section 3.2.1.

The solution is based on the construction of two independent light cones, as illustrated on the right panel of Fig. 5. The first light cone extends from the observer at z = 0 up to the redshift surface where the line of sight intersects the boundaries of the stack. This is referred to as the foreground light cone. The second light cone expands from the last redshift surface of the foreground light cone up to the maximum redshift, z_max, and has a different orientation compared to the first one to ensure that no box boundaries are hit out to z_max. This is referred to as the background light cone. The line of sight of each light cone is specified by the tuple of positions defined by the target object and the observer. For the foreground light cone, the target object is a specific selected halo in the simulation and the observer position is randomly generated on the XY plane. In the case of the background light cone, the observer is located underneath a randomly selected position of the last box of the stack. Each of the light cones are then assembled following the approach described in Section 3.2.1. Clearly the axes of symmetry are misaligned since they are built independently and point to different directions. Nevertheless, for the light-cone construction, the only relevant quantity is the relative angles of an object with respect to the axis of the light cone, i.e. δRA and δDec. These are independent of the direction of the light-cone axis. Hence, they can be used to align the two independent light cones, the foreground and background ones, to point to the same direction. We refer to this methodology as cut-and-paste. This method may affect the large-scale correlation function of AGN, although this is expected to be a second-order effect and does not modify our results and conclusions.

3.3 Implementation for this work: simulating a realistic set of observations

We use the implementation of universemachine (Behroozi et al. 2019) on the MDPL2 (Klypin et al. 2016) cosmological simulation with a side of 1 Gpc h⁻¹. We select a total of 12 universemachine boxes at different snapshot redshifts chosen to cover the redshift range z = 0 – 3 in steps of ∼1 Gyr. Mock AGNs are assigned to universemachine galaxies using the SARDs of either Georgakakis et al. (2017) or Aird et al. (2018). Our baseline simulations use as reference the observations of the cluster PLCKG266.6-27.3 with a mass of M|$_{500}^{\rm SZ}=8.5\times 10^{14}$| M_⊙ at a redshift of z = 0.97 (see Bartalucci et al. 2018). This is because PLCKG266.6-27.3 is the cluster in the sample of Koulouridis & Bartalucci (2019) that shows the highest excess of X-ray sources at a projected radial distance of 2.5 |${R}_{500}$|⁠. The mock Chandra/ACIS-I observations of PLCKG266.6-27.3 use massive haloes drawn from the universemachine box at a snapshot redshift of z = 0.94, i.e. similar to the redshift of the real cluster. There are 10 haloes in that box with M_500c>5 × 10¹⁴ M_⊙,² i.e. similar to the limiting mass of the Koulouridis & Bartalucci (2019) sample. The light cones are constructed to target the most massive halo in the simulation box with a mass of M_500c ∼ 7.5 × 10¹⁴ M_⊙ (universemachine identification number id = 7830644447). We study the impact of halo mass on our results and conclusions by also constructing light cones that pass through a second less massive halo (universemachine identification number id = 7793510527) with mass M_500c ∼ 5 × 10¹⁴ M_⊙. In the next section, we show that our analysis is not sensitive to the choice of the massive halo used to simulate light cones of clusters of galaxies.

We generate 100 lines of sights pointing to each of these two clusters with a field of view set to 20 arcmin diameter, which mimics the Chandra/ACIS-I observations. For each simulated observer we produce the projected radial distribution of mock X-ray selected AGN by splitting the field of view in eight concentric rings. The i-th ring is assigned an outer radius r_i = i · 0.5 R₅₀₀ from the cluster centre. Mock galaxies, and therefore, AGN of the light cone are assigned to a ring depending on their projected radial distance relative to the cluster centre. The application of observational selection effects onto the simulation requires the estimation of AGN fluxes. They are assigned to the X-ray AGN luminosities by assuming a power-law spectral shape with photon index Γ = 1.4, i.e. similar to that of the diffuse X-ray background at energies below about 10 keV (Akylas et al. 2012). Applying the corresponding sensitivity curve of each ring (see Section 2), a probability of detection is assigned to each source based on its flux. The total number of AGN per ring is calculated as the sum of probabilities of the AGN within the ring. The final product of this process are 100 AGN radial distributions that represent the 100 simulated lines of sight. The background is statistically subtracted as in observations (see Section 2). We calculate the expected number of background/foreground AGN within each ring by simulating field (i.e. off cluster) observations. We generate a field sample by constructing 100 light cones that point to 100 random locations in the box at z = 0.97, i.e. the same box as the simulated clusters. For simplicity we always locate the observer underneath the random target point. The same set of split overlap surfaces used for the clusters is also applied to the field observations. This is because we require the same redshift structure in both samples. We calculate the AGN distribution for this sample following the same steps described for the clusters. Each AGN is assigned to one ring using its projected radial distance with respect to the centre of the field. Detection probabilities are assigned to the AGN using the corresponding area curve. Finally, the expected field value is calculated as the average number of AGN per radial bin in the 100 simulations.

4 RESULTS

4.1 AGN halo occupation predictions

Before focusing on the radial distribution of AGN in massive clusters of galaxies at z ≈ 1, we present general predictions of our semi-empirical model (see Section 3.1) on the halo occupation distribution (HOD) of X-ray selected AGN. Such model predictions can be compared against current and future observational constraints to gain insights into the triggering mechanisms of accretion events onto SMBHs at different environments. We reiterate that our semi-empirical model is build upon the fundamental assumption that the clustering of AGN follows that of their host galaxies. The latter is included in the modelling of the galaxy–halo connection as implemented by universemacine. Any discrepancies between observed AGN HODs and the semi-empirical model predictions would question the assumption above, thereby pointing to environmental effects that modulate the incidence of AGN in haloes (e.g. Muñoz Rodríguez et al. 2023). We also remind the reader that the AGN–galaxy connection approach presented in Section 3.1 produces mock AGN catalogues that are consistent with the observed 2-point correlation function of different AGN samples that span a range of accretion luminosities and redshifts (Georgakakis et al. 2019). In that respect our semi-empirical model is consistent with the large-scale distribution of AGN in the Universe.

The AGN HOD, 〈N(L_X)|M〉, is defined as the mean number of AGN brighter than the luminosity L_X in dark matter haloes of given mass, M. Because of the different halo types (central or satellites) the HOD is usually expressed as a sum of two terms

where |$\langle N_{\rm cen}(\rm L_X)| M \rangle$|⁠, |$\langle N_{\rm sat}(\rm L_X)| M \rangle$| is the mean number of AGN brighter than L_X in parent haloes of mass M that are associated with central and satellite galaxies, respectively. f_A is a normalization factor that represents the fraction of active galaxies with respect to the full population.

Fig. 6 shows our HOD predictions for different X-ray luminosity cuts and redshifts. These are estimated by populating the universemachine boxes at the corresponding redshifts with AGN and then applying equation (2). At fixed luminosity threshold the HOD normalization increases towards higher redshift. This is the result of the strong increase of the AGN space density to redshift z ≈ 2–3. Also, the HOD normalization decreases towards higher luminosities as a result of the form of the AGN X-ray luminosity function. Fig. 6 further shows that both specific accretion rate models used to seed galaxies with AGN produce similar HOD results. However, the differences between the two models are stronger at the lowest redshift bin (z ∼ 0.25) and towards higher X-ray luminosities. These discrepancies are related to the modelling of the observed specific accretion rate distributions by Georgakakis et al. (2017) and Aird et al. (2018) as already discussed in Muñoz Rodríguez et al. (2023).

Fig. 6 also compares our semi-empirical model predictions with recent results on the HOD of X-ray selected AGN in the eROSITA eFEDS field (Comparat et al. 2023). This sample selects AGN in the redshift interval z = 0.05–0.55 (average of 0.34) and mean X-ray luminosity in the 0.5–2 keV band of |$\approx 10^{43}\rm erg\, s^{-1}$|⁠. Two clustering statistics, the 2-point correlation function and weak lensing, are applied to this sample to measure the AGN halo occupation distribution. We caution that the normalization of the AGN HOD, i.e. parameter f_A in equation (2), cannot be inferred from the observations (Allevato et al. 2021; Carraro et al. 2022). Instead this important quantity is determined post-processing based on knowledge of the AGN X-ray luminosity function and halo mass function at the redshifts of interest (e.g. Krumpe et al. 2023). For the comparison we fixed f_A = 0.01, which correspond to the duty cycle of central galaxies derived from the specific accretion rate distributions at similar redshift and luminosity threshold to those used by Comparat et al. (2023). Although the uncertainties of the observations are large, there is evidence that the best-fitting AGN HOD increases with increasing halo mass slower (i.e. flatter slope) than the model predictions. This is consistent with claims for suppression of AGN activity towards massive haloes, i.e. clusters of galaxies, at z ∼ 0.25 compared to the less dense regions of the cosmic web, e.g. field (e.g. Haines et al. 2012; Mishra & Dai 2020). This is also in line with the arguments presented in Muñoz Rodríguez et al. (2023) based on the forward modelling of the observed fraction of AGN in massive clusters of galaxies.

4.2 Observed overdensity of AGN

Next we compare the projected radial distribution of X-ray selected AGN in the cluster PLCKG266.6-27.3 (see Fig. 1) with the predictions of the semi-empirical model described in the previous sections. Fig. 7 shows this comparison for two versions of the model based on either the Aird et al. (2018) or the Georgakakis et al. (2017) specific accretion rate distributions. The MDPL2 halo selected to represent the cluster PLCKG266.6-27.3 has a catalogued identification number id = 7 830 644 447 in universemachine and a halo mass of M_500c ∼ 7.5 × 10¹⁴ M_⊙. At fixed R₅₀₀ radial bin, Fig. 7 plots the mean of the model predictions and the corresponding 68 per cent confidence intervals. These quantities are determined from the radial distributions of individual fiducial observers. The scatter around the mean (68 per cent confidence interval) therefore provides an estimate of the sample variance, i.e. the fact that different observers see different structures along their corresponding lines of sight to the cluster.

The simulations predict, on average, a flat radial distribution independent of the adopted specific accretion rate model used to seed galaxies with AGN. This is an expected behaviour of the model, which assumes that the incidence of AGN in galaxies (i.e. the probability of triggering an accretion event onto a SMBH) is independent of environment. As a result, there is no special physical scale in the model at which an overdensity of AGN should be expected. A striking feature in Fig. 7 is the large scatter around the mean at fixed R₅₀₀ radial bin. In that respect, it is interesting that within the 1.5σ sample variance uncertainty the predictions of the models are in agreement with the observations. We reiterate that the origin of this scatter is the diversity of projected structures along the line of sight of different observers. It is therefore interesting to explore whether individual simulated observers (i.e. individual light-cone realizations) see X-ray AGN radial distributions with features similar to the observed ones, i.e. excess counts.

Fig. 8 shows the radial distributions for each of the 100 fiducial observers. Eyeballing each of these realizations would identify a few that show an X-ray AGN overdensity at the radial distance of 2.5 R₅₀₀ relative to the neighbouring bins. This approach however, is subjective and therefore we define a set of quantitative criteria to select light cones with excess counts. The adopted conditions that should be simultaneously met are

where N_AGN(r_i) is the number of AGN at the ring i, i indicates the ring of the overdensity (i.e. r = 2 – 2.5 R₅₀₀), σ_i is the scatter in the correspondent radial bin, and i ± 1 indicate the previous and subsequent ring (i.e. r = 1.5 – 2 R₅₀₀ and r = 2.5 – 3 R₅₀₀, respectively). We reiterate that this criterion is empirically motivated, i.e. it is tuned to broadly select simulated radial distributions similar to the observations. Therefore, visual inspection of Fig. 8 may reveal either simulated radial distributions that fulfil the criteria but show marginally significant peaks at r = 2 – 2.5 R₅₀₀ (e.g. lc = 14, 38, or 47) or, conversely, realizations that show excess counts at that ring but are not picked by the criteria (e.g. lc = 3, 6, or 44). We acknowledge these issues, which on the other hand, emphasize the necessity of having a quantitative, objective, and reproducible approach for selecting simulated projected radial distributions. Hence, equation (3) provides a basis for the quantitative assessment of the frequency of AGN overdensities in their projected radial distribution. Fig. 8 highlights the realizations that fulfil the above criteria. It demonstrates that ∼20 per cent of the observers reproduce similar peaks as in the observations. This frequency is only mildly sensitive to the adopted criteria. We therefore conclude that sample variance needs to be taken into account when interpreting the radial distribution of AGN in massive clusters of galaxies. A non-negligible fraction of our simulation realizations can reproduce the most extreme cluster, in terms of excess AGN, of the sample of Koulouridis & Bartalucci (2019).

For the simulated observations in Fig. 8 that reproduce an excess of X-ray counts at |$2.5\, R_{500}$|⁠, we further explore the redshift distribution of mock AGN. This is to investigate if the excess of X-ray sources is associated with the cluster. In Fig. 9 we show seven examples of the redshift distribution of mock AGN in the radial distance bin of |$2.5\, R_{500}$|⁠. These are selected from light-cone configurations in Fig. 8 that reproduce an overdensity of mock AGN at |$2.5\, R_{500}$|⁠. Most of these realizations show a redshift distribution where the peak is generated by objects in the foreground and/or the background of the cluster. There are also realizations where mock AGN that produce the overdensity are at redshifts similar to the cluster. We reiterate, however, that even in this case this is a projection effect because of the zero-order assumption of the model. Nevertheless, the contribution of these cases is marginal and most of the redshift distributions are dominated by foreground and/or background sources.

Next, we explore the incidence of excess projected X-ray counts in other R₅₀₀ rings around the simulated clusters. We adopt the same set of conditions defined above to identify in a quantitative manner excess counts. The only deviation is that the main ring within which overdensities are searched for varies between r = 1.5–3.5 R₅₀₀. The model predicts an occurrence of about 10–20 per cent of an overdensity for different cluster centric distances. For distances smaller than r = 1.5 R₅₀₀ or bigger than r = 3.5 R₅₀₀, the sensitivity of the observation drops because of the extended emission of the cluster and the increasing off-axis incidence angles, respectively. Hence, it is difficult to make a clear comparison at these radii.

All the results above correspond to simulations of a single massive halo (M_500c ∼ 7.5 × 10¹⁴ M_⊙) and the implementation of a single sensitivity map, the one that corresponds to the Chandra observations of PLCKG266.6-27.3 with a total on source exposure of ≈200 ks. Next we explore the impact of different X-ray depths in the result and conclusions. For this purpose, we repeat the analysis using the same halo in the simulations to construct light cones but, applying a different sensitivity map to construct mock X-ray observations. The new map corresponds to the shallower Chandra observations of the cluster SPTCLJ2146-4633 with a total on source exposure of ≈70 ks. The main effect is that the number of detected AGN and the overall scatter decreases. This is because less luminous AGN are harder to detect in the case of shallower X-ray observations. Nevertheless, since the total number of AGN also decreases the fraction of mock observers that see an excess of projected X-ray counts at |$r=2.5\, R_{500}$| based on the conditions presented earlier (see equation 3) is similar to our baseline results using the most sensitive observation, i.e. about 10 per cent (see also upper panels of Fig. 10). The effects of cluster mass onto the radial distribution are also investigated by repeating the same exercise for a different less massive halo in the N-body simulation with M_500c ∼ 5 × 10¹⁴ M_⊙ (universemachine id = 7793510527). The corresponding radial distributions for the different sensitivity maps (i.e. the ones of the observed clusters PLCKG266.6-27.3 and SPTCLJ2146-4633) are shown in the lower panels of Fig. 10. In both cases we find a flat mean distribution with a large scatter around it which mimics our baseline result. This is an expected feature of the model since its zero-order assumption is that the AGN activation is independent of the halo mass.

5 DISCUSSION

5.1 Radial distribution of AGN

The overarching question of this work relates to the role of the environment in modulating accretion events onto the SMBHs at the nuclear regions of galaxies. We approach this problem by investigating the X-ray AGN projected radial distribution in the vicinity of massive clusters of galaxies. These structures represent the densest regions in the Universe, where environmental effects and processes are expected to reach their maximum impact (e.g. starvation, strangulation, or ram-pressure, see Gunn & Gott 1972; Larson, Tinsley & Caldwell 1980; Moore et al. 1996). It is now well established that the number of AGNs, in clusters of galaxies increases with redshift (e.g. Martini et al. 2009; Martini et al. 2013). This trend mirrors the evolution of the overall AGN field population (e.g. Ueda et al. 2014; Aird et al. 2015) and perhaps proceeds even faster (e.g. Ehlert et al. 2014; Bufanda et al. 2017; Hashiguchi et al. 2023; Toba et al. 2024), thereby suggesting that dense environments at high redshift promote accretion events onto SMBHs (e.g. Lehmer et al. 2009; Digby-North et al. 2010). It has been proposed that the incidence of AGN in massive clusters is related to an infalling population of galaxies whose black holes become active as they enter the dense cluster environment (e.g. Haines et al. 2012; Pimbblet et al. 2013; Rihtaršič et al. 2024).

In this work, we test this scenario by modelling the observed projected radial distribution of X-ray selected AGN in massive clusters at z ≈ 1 presented by Koulouridis & Bartalucci (2019). That cluster sample is advantageous because the individual cluster properties (mass and radius) are accurately determined using a sophisticated method that combines information from both XMM–Newton and Chandra observations (see Bartalucci et al. 2017, 2018). The large effective area of the former allows the characterization of faint structures, while the spatial resolution of the latter enables modelling the central regions of the clusters. This approach leads to an accurate characterization of the density profile of the clusters out to R₅₀₀. For this sample it is therefore possible to build robust radial distributions of X-ray selected AGN as function of distance normalized to |${R}_{500}$| and explore evidence for a statistically significant excess of counts at the radius |$2.5\, {R}_{500}$|⁠.

The semi-empirical modelling developed in this work emphasizes the role of sample variance in the interpretation of the observed projected AGN radial distributions. We produce mock AGN catalogues under the explicit assumption that accretion events on the SMBHs are triggered with the same probability in the different environments. Then we use these mock AGN and galaxy catalogues to simulate realistic observations of clusters that include the same selection effects as the observations of Koulouridis & Bartalucci (2019). We study the impact of projection effects by simulating 100 observations of the same cluster in the simulation with randomly selected lines of sight (see Section 3.3). A striking results from our analysis is the flatness of the simulated average projected radial distribution (see Fig. 7), which at first instance appears inconsistent with the observations of Koulouridis & Bartalucci (2019). At the same time however, there is substantial scatter around the mean of this distribution as a result of sample variance, i.e. background/foreground structures along the line of sight projecting into the field of view (see Fig. 9). Given this scatter the significance of the excess counts at the radial ring |$2.5\, {R}_{500}$| in Fig. 7 is significant only at the |$1-2\, \sigma$| level. We nevertheless, take a further step and calculate the probability of finding overdensities similar to the observed ones. This analysis also demonstrates the importance of stochasticity in producing excess X-ray AGN counts in the radial distribution of counts that have no physical origin. The model reproduces radial distribution overdensities at |$2.5\, {R}_{500}$| similar to those found by Koulouridis & Bartalucci (2019) in up to 20 per cent of the simulated light cones (see Fig. 8). This fraction should be compared with the rate of 40 ± 20 per cent (two out of a total of five clusters, we assume binomial statistics for the estimation of the uncertainty in this fraction) in the sample of Koulouridis & Bartalucci (2019) that show a statistically significant excess of counts. These results also have implications for other studies in the literature that find evidence for excess AGN counts in the projected radial distribution of AGN beyond the virial radius (Johnson et al. 2003; Ruderman & Ebeling 2005; Fassbender et al.2012).

We caution that our simulations cannot reject the possibility of a physical interpretation of the excess counts at |$2.5\, {R}_{500}$| found by Koulouridis & Bartalucci (2019). For example Rihtaršič et al. (2024) suggest, based on hydrodynamical simulations, a higher fraction of AGN (⁠|$L_{\rm X}(0.5-10\, \rm {keV})\gtrsim 10^{42}\, \rm {erg}\,\rm {s^{-1}}$|⁠) among massive galaxies (M_⋆ ≳ 10¹¹ M_⊙) in groups or clusters at radial distances of about 3 R₅₀₀. They interpret this result as evidence for a preferential activation of black holes in the group/cluster outskirts, but caution that projection effects because of substructure (e.g. infalling groups) may swamp this signal. Addressing the origin of this excess, physical or stochastic, requires spectroscopic information, which would allow the robust identification of AGN cluster members and separate them from foreground/background interlopers. Increasing the cluster sample will allow a better understanding of the physics at play. This is because different studies show that the dynamical state (i.e. relaxed versus non-relaxed) of the cluster could have an impact on the AGN activity (see Kocevski et al. 2009; van Breukelen et al. 2009; Stroe & Sobral 2021) and the two clusters which show the overdensity in Koulouridis & Bartalucci (2019) are in different relaxation states, i.e. one of them is virialized while the other is not (see Bartalucci et al. 2017).

5.2 Exploring a higher incidence of AGN among the infalling mock galaxy population

Next, we test the hypothesis of an infalling population as the origin of the excess counts in the radial distribution of AGN at about |$2\, {R}_{500}$| in Fig. 7. Our approach is to tune our semi-empirical model by associating a higher incidence of AGN among infalling galaxies. This requires (i) a criterion for isolating galaxies that enter for the first time the cluster from the cosmic web and (ii) a new specific accretion rate distribution model that is applied to these galaxies and has the property of producing a higher incidence of AGN at fixed X-ray luminosity threshold.

Ideally, an infall population would be defined by following the orbits of the dark matter particles that make up haloes. However, semi-empirical models, like universemachine, are build upon halo catalogues and therefore information about the formation/merging history of individual haloes is not readily available. Instead, we decide to adopt the alternative but widely used approach of the phase-space diagram to find infalling haloes. For a given massive cluster halo in universemachine it is possible to estimate the relative velocities (v_3D) and relative radial distances (⁠|$\rm {r}_{\rm 3D}$|⁠) to other haloes in the simulation (parent or satellites). The phase-space diagram of the cluster under consideration is then defined by the parameters v_3D/σ_cl (σ_cl is the velocity dispersion of the main cluster halo) and r|$_{\rm 3D}/{R}_{500}$| (R₅₀₀ refers to the main cluster halo). Haloes with different infall histories populate distinct regions of the phase-space plane. This is because the ratio between r|$_{\rm 3D}/{R}_{500}$| and v_3D/σ_cl is a proxy of the infall time of a main cluster halo (see e.g. Noble et al. 2013, 2016; Rhee et al. 2017; Kim et al. 2023). In the notation above the |$3\rm {D}$| index refers to 3D quantities estimated from the spatial distribution of haloes in the universemachine simulation box. The adopted |$3\rm {D}$| phase-space diagram is independent of projection effects that are inevitable when constructing light-cone realizations assuming different observer positions (see Section 3.2). Following commonly used criteria we define infalling haloes/galaxies as those or that simultaneously satisfy the following conditions

where |$v_{esc,\rm {NFW}}$| corresponds to the escape velocity (e.g. Rhee et al. 2017) of a halo assuming an Navarro–Frenk–White profile (NFW; Navarro, Frenk & White 1996) and |$\rm {M}_{halo,peak}$| is the maximum historical mass of the halo. Fig. 11 shows the phase space diagram for the cluster with dark matter halo id = 7830644447, i.e. the same massive halo used to construct light cones and simulated radial distributions (see Section 3.3). The sample of infalling galaxies based on the conditions above is indicated with the red circles in Fig. 11.

The next step is to adopt a new specific accretion rate distribution model, which when applied to the infalling galaxies above yields a higher fraction of AGN. In Muñoz Rodríguez et al. (2023) we showed that the observed fraction of X-ray selected AGN relative to galaxies in massive clusters of galaxies at z ≈ 1 is much higher than that predicted by our baseline semi-empirical model that uses either the Georgakakis et al. (2017) or the Aird et al. (2018) SARs. Instead, Muñoz Rodríguez et al. (2023) proposed that a log-normal SAR model with mean specific accretion rate log λ_SAR = − 1.25 and scatter σ = 0.1 applied to galaxies with stellar masses M_* > 10^10.7 M_⊙ can reconcile the tension with the observed fraction of X-ray selected AGN in massive clusters at z ≈ 1. We therefore choose to use the same SAR model in our analysis and apply it only to the infalling galaxies (red circles) of Fig. 11. The impact on the AGN radial distribution of the modified SAR for the infalling galaxies is shown in Fig. 12. Relative to our baseline model the mean expected number of X-ray selected AGN slightly increases for the radial ring |$2-2.5\, {R}_{500}$|⁠, i.e. the one where excess counts where observed by Koulouridis & Bartalucci (2019). However, the same effect is seen at smaller radii, |$0.5-2\, {R}_{500}$|⁠. This is because the infall population is evenly distributed between |${r}_{\rm 3D}=0.5-3\, {R}_{500}$| as demonstrated by the top panel of Fig. 11. In any case, the increase at the ring |$2-2.5\, {R}_{500}$| is modest and is associated with substantial scatter. We apply the criteria of equation (3) to identify in an objective manner excess counts in the ring |$2-2.5\, {R}_{500}$| among the light-cone realizations with the modified SAR. We find that 20 per cent of the light cones show radial distributions that resemble the observations. This rate is the same as with the baseline semi-empirical model predictions presented in Fig. 8. We conclude that the approach outlined above for increasing the incidence of AGN among infalling galaxies in massive clusters has a moderate impact on the observed radial distribution of AGN and cannot fundamentally modify the predictions of our baseline semi-empirical model.

6 CONCLUSIONS

In this paper we develop a flexible semi-empirical model of AGN and galaxies in a cosmological volume to interpret observations of the radial distribution of AGN in massive clusters of galaxies at z ≈ 1 (Koulouridis & Bartalucci 2019) and test claims for an efficient activation of SMBHs in the outskirts of galaxy clusters. The explicit assumption of the model is that the AGN triggering is independent of environment (or halo mass). This allows us to test the hypothesis that the excess counts of X-ray selected AGN observed at a radius of about |$2-2.5\, {R}_{500}$| in massive clusters of galaxies at z ≈ 1 (Koulouridis & Bartalucci 2019) are not physical but instead driven by projection effects.

We select haloes at z ≈ 1 in the simulations with masses similar to the clusters of Koulouridis & Bartalucci (2019) and generate mock observations through different sightlines to test the impact of sample variance to the inferred mock AGN radial distribution. A key step of this process is the generation of light cones which allows us to implement the selection effects of the real observations to the mocks (e.g. field-of-view, variations of the flux limit at different radial distances from the cluster centre). The main results of the paper are:

• We demonstrate that our semi-empirical model predicts HODs for X-ray selected AGN in broad agreement with the latest observational constraints of Comparat et al. (2023) at z ∼ 0.2. The normalization of our HODs decreases towards lower redshift and brighter luminosities, mirroring the evolution of the X-ray AGN population with redshift and the form of the X-ray luminosity function.

• There is evidence for a possible tension between observations and model predictions on the HOD slope of satellite AGN. The observations favour flatter slopes compared to the semi-empirical model. Although the observational uncertainties are large, this discrepancy may point to the suppression of X-ray AGN in satellites galaxies of massive cluster of galaxies at z ≈ 0.25.

• Turning to the projected radial distribution of X-ray selected AGN in the vicinity of massive clusters at z ≈ 1, our model predicts on average a flat radial distribution. This is a direct consequence of the main assumption of the model construction that the AGN triggering is independent of the environment (Fig. 7).

• Our analysis emphasizes the importance of sample variance that manifests as scatter around the mean of the projected radial distributions predicted by the model. As a result in a non-negligible number of cases excess counts at radial distances of 2–2.5 R₅₀₀ are predicted by the model. Up to 20 per cent of the realizations show amplitudes similar to the observations of Koulouridis & Bartalucci (2019) for massive clusters of galaxies at z ≈ 1 (see Fig. 8). This incidence rate is lower but still consistent within the errors with the observed fraction of clusters in the Koulouridis & Bartalucci (2019) work with excess counts in their outskirts, 40 ± 20 per cent. In our model, however, these overdensities in the projected radial distribution are not physical but stochastic and dominated by interlopers (Fig. 9).

• Fine tuning our model to favour a higher incidence of mock AGN among galaxies in the infall region of massive haloes has little impact to the predicted projected radial distributions (see Fig. 12). This further emphasizes the significance of sample variance in interpreting projected AGN radial distributions.

Funding

This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 860744. AG and IMR acknowledge support from the Hellenic Foundation for Research and Innovation (HFRI) project ‘4MOVE-U’ grant agreement 2688, which is part of the programme ‘2nd Call for HFRI Research Projects to support Faculty Members and Researchers’. AL is partly supported by: ‘Data Science methods for MultiMessenger Astrophysics & Multi-Survey Cosmology’, Programmazione triennale 2021/2023 (DM n.2503 dd. 9 December 2019), Programma Congiunto Scuole; Fondazione ICSC, Spoke 3 Astrophysics and Cosmos Observations Project ID CN-00000013.

ACKNOWLEDGEMENTS

The authors would like to thank the anonymous referee for their comments to the paper. 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), NASA’s Astrophysics Data System (ADS) and the arXiv preprint server.

DATA AVAILABILITY

The data set of the simulated light cones used in the analysis are available in Zenodo⁸ and the code to reproduce the results of the paper is available in a Github repository.⁹ The light-cone generation code is available upon request and will be public on Github in the near future.

Footnotes

REFERENCES