Scrutinising evidence for the triggering of Active Galactic Nuclei in the outskirts of massive galaxy clusters at $z\approx1$

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; AGN). 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 halos at that low redshift. We then use this model to interpret the projected radial distribution of X-ray sources around high redshift ($z\approx1$) massive ($>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% 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 emphasises the importance of projection effects and shows that current observations of $z\approx1$ clusters remain inconclusive on the activation of SMBHs during infall.


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 realisation that nearly every spheroidal galaxy hosts at its nuclear regions a black hole with a mass that may exceed 10 9 M ⊙ 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 sug-★ E-mail: ivan.rodriguez@noa.grgest 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 (AGN).
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).Disk instabilities (e.g.Hopkins & Hernquist 2006;Gatti et al. 2016) and galactic bars (Cis-ternas et al. 2015) are efficient in removing angular momentum from the interstellar gas thereby driving it toward 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 et al. 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 et al. 2005;Koulouridis et al. 2006Koulouridis et al. , 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.Observations of the fraction of AGN in low redshift ( < ∼ 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;Sabater et al. 2013;Martini et al. 2013;Koulouridis et al. 2014;Mishra & Dai 2020).Nonetheless, this trend seems to level off at intermediate redshift ( ≈ 0.7, Eastman et al. 2007;Galametz et al. 2009;Martini et al. 2009Martini et al. , 2013;;Ehlert et al. 2014;Bufanda et al. 2017) and perhaps reverse at  > ∼ 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  ≈ 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.There is indeed evidence that the fraction of AGN relative to galaxies is decreasing toward 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 et al. 2017;Mishra & Dai 2020;Stroe & Sobral 2021;Koulouridis et al. 2024).Additionally, there are claims that the projected counts of AGN show an excess outside the virial radius of clusters (Johnson et al. 2003;Ruderman & Ebeling 2005;Fassbender et al. 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 et al. 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  500 from the cluster centre, where  500 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  500 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 halos 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 Ω  = 0.307, Ω Λ = 0.693, ℎ = 0.678 consistent with the Planck results (Planck Collaboration et al. 2016).

OBSERVATIONS
This work uses Chandra X-ray observations of massive clusters of galaxies at  ≈ 1 presented by Koulouridis & Bartalucci (2019).Their sample is selected using the Sunyaev-Zeldovich effect (SZ) and it is composed by the 5 most massive clusters (  500 ≳ 5 × 10 14 M ⊙ ) in the South Pole Telescope and Planck catalogues at that redshift (see Bleem et al. 2015;Planck Collaboration et al. 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  500 (0.7 − 1 Mpc) and M 500 (mass range 3 − 8 × 10 14  ⊙ ).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 ≈ 2.5  500 .In this paper, we build a semiempirical model of the radial distribution of AGN in massive dark matter halos 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 datasets 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.0keV (full), 0.5-2.0keV (soft), 2.0-7.0 keV (hard) and 4.0-7.0keV (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 −4 .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 −6 .X-ray fluxes are determined assuming a power-law X-ray spectrum with photon index Γ = 1.4 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 one-dimensional 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 Xray 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 •  500 .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  − log  distribution.For a given cluster and  500 radial bin, , we first determine the full-band sensitivity curve of the ring with inner and outer radius of /2 •  500 and ( + 1)/2 •  500 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 Figure 1.The two selected clusters are the ones that show the highest overdensity at 2.5  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 and SPTCLJ2146-4633 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 dataset, ≈ 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.

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 et al. 2018;Behroozi et al. 2019;Grylls et al. 2019) and/or AGN (e.g.Comparat et al. 2019Comparat et al. , 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 Dark matter halos are populated with galaxies using abundance matching techniques (panel labelled stellar mass).Supermassive black-hole accretion events are painted of top of these galaxies using empirical relations that describe the probability of a galaxy hosting an AGN (panel labelled AGN).See Section 3.1 for details on the mock catalogue construction.The right panel on the upper branch represents a light cone with a field of view 20 arcmin pointing to a massive halo (i.e.cluster of galaxies) in the mock catalogue as described in Section 3.2.The lower branch shows the derivation of the X-ray selection effects.The starting point are X-ray observations (left panel) of cluster of galaxies presented in Section 2. These are processed to obtain the X-ray sensitivity maps (middle panel) explained in Section 2. The anulii in these panels have radii that are multiples of the quantity 0.5 •  500 .These are the anulii used in the observations (see Figure 1) to study the radial distribution of AGN in clusters.The different colours represent different rings.The right panel shows the area curve for the different rings which describe the probability of detecting a source at different radial distances from the centre of the cluster.The colour of each line corresponds to the anulii shown in the sensitivity map panel.Both branches merge in the panel labelled "Simulated observation", which represents a simulated field of view that includes the selection effects derived from observations.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/ℎ, a mass resolution of 1.5 × 10 9 M ⊙ /ℎ and 3 840 3 (∼57•10 9 ) particles.Dark matter halos 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 halos by parameterising 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 halos, 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; Bongiorno et al. 2012;Aird et al. 2012;Bongiorno et al. 2016;Georgakakis et al. 2017;Aird et al. 2018).This quantity describes the probability of a galaxy hosting an accretion event onto its supermassive black hole with specific accretion rate (SAR)    =   / ★ , where   is the AGN luminosity (in this case at X-rays) and  ★ 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 Figure 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 the present 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.

Light-cones
In this work, we explore the projected radial distribution of X-rayselected 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 halos whose sight lines probe different paths through the cosmic web.This is demonstrated in Figure 3 which shows two different sight lines to a particular halo (left panel).The corresponding projected structures along these sight lines 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 sight lines to a particular halo (Section 3.2.2).

General light cone construction
Extragalactic surveys are typically characterised 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 1 .Producing mocks over a wide redshift interval requires that simulated boxes are used as building blocks to construct a 3-dimensional (3D) pavement.The stack of boxes can be extended from  = 0 to an arbitrary maximum redshift (  ) 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.
1 For example, the centre of a box with a length-size of 1 Gpc/ℎ at  = 1 corresponds to a comoving distance of  , ∼ 2300 Mpc.The bottom and top faces lie at  , ∼ 1800 Mpc and  ,  ∼ 2800 Mpc respectively.These distances correspond to the redshift range  ∼ 0.72 − 1.3.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).
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  /2, L  /2], where L  is the length-side of the simulated volume.The hypothetical observer is located onto the XY plane at  = 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 (    ) with respect to the position of the observer.This is achieved by offsetting the box along the Z-axis by Δ defined as where   ( =     ) is the comoving distance that corresponds to     .The coordinates in the equation above are defined as x = x  x  , y = y   -y  and z = z   , where  are the coordinates of the observer and  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   ( =     ) 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     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   (  ) −   ( +1 ) <   , 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 Figure 4, which shows the stack of two boxes.The intersecting volume contains the individual halos 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 ().It determines which objects are adopted from each box.Above the surface, only halos 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 Figure 4, where the lower curved line represents the splitoverlap-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 sight line 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 sight line is calculated.This angle can be decomposed into a right-ascension and declination on the unit sphere.The redshifts associated with the individual halos 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 Figure 4, where only objects within the limits of the field of view are included.

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 Figure 5.The sight line 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 Figure 5.The first light cone extends from the observer at  = 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,   , and has a different orientation compared to the first one to ensure that no box boundaries are hit out to   .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.   and .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.

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/ℎ.We select a total of 12 UniverseMachine boxes at different snapshot redshifts chosen to cover the redshift range  = 0 − 3 in steps of ∼ 1 Gyr.Mock AGN 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   500 = 8.5×10 14 M ⊙ at a redshift of  = 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  500 .The mock Chandra/ACIS-I observations of PLCKG266.6-27.3use massive halos drawn from the UniverseMachine box at a snapshot redshift of  = 0.94, i.e., similar to the redshift of the real cluster.There are 10 A sketch that illustrates the problem of using arbitrary lines of sight when constructing light cones that are forced to intercept a particular halo.The panels show a stack of 3 boxes at different snapshot redshifts.Each of the boxes is coloured differently (blue, green or red) for illustration purposes.The blue and red shadings correspond to the lowest and highest redshifts of the stack.The position of the observer is shown at the bottom of each stack of boxes with the eyeball graphic.The left panel shows an example of a tilted sight line that is forced to contain the position of a halo marked with the star symbol in the green box.This sight line hits the boundaries of the stack of boxes before reaching the expected maximum redshift   (see Section 3.2.1).Such a light cone is clearly incomplete.The right panel visualises a solution to this issue that is based on the construction of two independent light cones.The first one encompasses the region between the observer and the last redshift slice where the light cone is complete.We refer to this component as the foreground light cone.A second independent light cone is then constructed that extends from the previous complete redshift to   .We refer to this component as background light cone.Finally, both light cones are aligned by matching the lines of sight.This is indicated in the far right panel.The arrow shows the rotation that needs to be applied to the background light cone to align it to the foreground one.
halos in that box with M 500 >5 × 10 14 M ⊙ 2 , 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 500 ∼ 7.5×10 14 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 500 ∼ 5 × 10 14 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 8 concentric rings.The -th ring is assigned an outer radius   =  • 0.5  500 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.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 2 UniverseMachine provides only virial halos masses.We convert these values to 500 critical, using a mean halo concentration  = 0.7 (Ludlow et al. 2014).
to 100 random locations in the box at  = 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.

AGN Halo occupation predictions
Before focusing on the radial distribution of AGN in massive clusters of galaxies at  ≈ 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 semiempirical 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 halos (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) where ⟨  (  )|⟩, ⟨  (  )|⟩ is the mean number of AGN brighter than   in parent halos of mass  that are associated with central and satellite galaxies respectively.  is a normalisation factor that represents the fraction of active galaxies with respect to the full population.
Figure 6 shows our HOD predictions for different X-ray luminosity cuts and redshifts.These are estimated by populating the Uni-verseMachine boxes at the corresponding redshifts with AGN and then applying Equation 2. At fixed luminosity threshold the HOD normalisation increases towards higher redshift.This is the result of the strong increase of the AGN space density to redshift  ≈ 2 − 3. Also, the HOD normalisation decreases toward higher luminosities as a result of the form of the AGN X-ray luminosity function.Figure 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 ( ∼ 0.25) and toward 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).
Figure 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  = 0.05 − 0.55 (average of 0.34) and mean X-ray luminosity in the 0.5-2 keV band of ≈ 10 43 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 normalisation of the AGN HOD, i.e., parameter   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   = 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, the best-fit AGN HOD increases with increasing halo mass slower (i.e.flatter slope) than the model predictions.This suggests a suppression of AGN activity toward massive halos, i.e., cluster of galaxies, at  ≈ 0.25 compared to the semi-empirical model expectation.This is 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.In that study it is found that the same semi-empirical model presented in Section 3.1 overpredicts the incidence of AGN among galaxies in low redshift clusters.

Observed overdensity of AGN
Next we compare the projected radial distribution of X-ray selected AGN in the cluster PLCKG266.6-27.3(see Figure 1) with the predictions of the semi-empirical model described in the previous sections.
Figure 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 = 7830644447 in UniverseMachine and a halo mass of M 500 ∼ 7.5 × 10 14 M ⊙ .At fixed  500 radial bin, Figure 7 plots the mean of the model predictions and the corresponding 68% confidence intervals.These quantities are determined from the radial distributions of individual fiducial observers.The scatter around the mean (68% 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 Cluster id=7830644447  2018) SARDs, respectively.The average at each radial bin is estimated from the 100 light cone realisations described in Section 3.3, which point to the massive halo ( 500 ≈ 7.5 × 10 14 M ⊙ /ℎ) with id=7830644447 in UniverseMachine.The light green and light orange shaded regions correspond to the 68% confidence intervals around the mean value at each radial bin.This scatter represents the (cosmic) variance among the 100 light cone realisations.The vertical lines represent the correspondent  200, (dashed) and   (dashed-dotted) normalised to the  500 of the cluster.
striking feature in Figure 7 is the large scatter around the mean at fixed  500 radial bin.In that respect, it is interesting that within the 1 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 realisations) see X-ray AGN radial distributions with features similar to the observed ones, i.e. excess counts.
Figure 8 shows the radial distributions for each of the 100 fiducial observers.Eyeballing each of these realisations would identify a few that show an X-ray AGN overdensity at the radial distance of 2.5  500 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    (  ) is the number of AGN at the ring ,  indicates the ring of the overdensity (i.e. = 2 − 2.5  500 ),   is the scatter in the correspondent radial bin, and  ± 1 indicate the previous and subsequent ring (i.e.,  = 1.5 − 2  500 and  = 2.5 − 3  500 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 Figure 9 may reveal either simulated radial distributions that fulfil the criteria but show marginally significant peaks at r=2-2.5  500 (e.g., lc=14, 38 or 47) or, conversely, realisations 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, emphasise 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.Figure 8 highlights the realisations that fulfil the above criteria.It demonstrates that ∼20% 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 realisations can reproduce the most extreme cluster, in terms of excess AGN, of the sample of Koulouridis & Bartalucci (2019).
For the simulated observations in Figure 8 that reproduce an excess of X-ray counts at 2.5  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 Figure 9 we show seven examples of the redshift distribution of mock AGN in the radial distance bin of 2.5  500 .These are selected from light cone configurations in Figure 8 that reproduce an overdensity of mock AGN at 2.5  500 .Most of these realisations show a redshift distribution where the peak is generated by objects in the foreground and/or the background of the cluster.There are also realisations 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  500 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  = 1.5 − 3.5  500 .The model predicts an occurrence of about 10-20% of an overdensity for different cluster centric distances.For distances smaller than  = 1.5  500 or bigger than  = 3.5  500 , 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 500 ∼ 7.5 × 10 14 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  = 2.5  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% (see also upper panels of Figure 10).The effects of cluster mass  2017) to seed galaxies with AGN.The observed projected radial distribution of X-ray selected AGN of the cluster PLCKG266.6-27.3 is shown with the grey/black squares connected with the solid grey/black lines (see Figure 1 and Section 2).The light-cone realisations that reproduce the observed peak at the distance of 2.5  500 are highlighted by (i) making the observational data points and connecting lines black, (ii) using bold green characters for the light-cone incremental number at the top of the corresponding panel and (iii) change in the background colour from white to grey.
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 500 ∼ 5 × 10 14 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 Figure 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.

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 Moore et al. 1996;Larson et al. 1980;Gunn & Gott 1972).It is now well established that the number of AGN in clusters of galaxies increases with redshift (e.g.Martini et al. 2009Martini 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).
In this work we test this scenario by modelling the observed projected radial distribution of X-ray selected AGN in massive clusters at  ≈ 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. 2017Bartalucci et al. , 2018)).The large effective area of the former allows the characterisation of faint structures, while the spatial resolution of the latter enables modelling the central regions of the clusters.This approach leads to an accurate characterisation of the density profile of the clusters out to  500 .For this sample it is therefore possible to build robust radial distributions of X-ray selected AGN as function distance normalised to  500 and explore evidence for a statistically significant excess of counts at the radius 2.5  500 .
The semi-empirical modelling developed in this work emphasises 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 Figure 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  500 in Figure 7 is significant only at the 1 − 2  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  500 similar to those found by Koulouridis & Bartalucci (2019) in up to 20% of the simulated light cones (see Figure 8).This fraction should be compared with the rate of 40±20% (2 out of a total of 5 clusters) 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  500 found by Koulouridis & Bartalucci (2019).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 vs. 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 virialised while the other is not (see Bartalucci et al. 2017).

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  500 in Figure 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 halos.However, semi-empirical models, like UniverseMachine, are build upon halo The upper row of panels is for the halo with id=7830644447 and mass  500 ≈ 7.5 × 10 14 M ⊙ /ℎ (same as in Figures 7, 8).The lower row of panels is for the halo with id=7793510527 and mass  500 ≈ 5 × 10 14 M ⊙ /ℎ (see Figure 1).
catalogues and therefore information about the formation/merging history of individual halos is not readily available.Instead, we decide to adopt the alternative but widely used approach of the phase-space diagram to find infalling halos.For a given massive cluster halo in UniverseMachine it is possible to estimate the relative velocities ( 3 ) and relative radial distances ( 3 ) to other halos in the simulation (parent or satellites).The phase-space diagram of the cluster under consideration is then defined by the parameters  3 /  (  is the velocity dispersion of the main cluster halo) and  3 / 500 ( 500 refers to the main cluster halo).Halos with different infall histories populate distinct regions of the phase-space plane.This is because the ratio between  3 / 500 and  3 /  is a proxy of the infall time of a main cluster halo (see e.g., Noble et al. 2013Noble et al. , 2016;;Rhee et al. 2017;Kim et al. 2023).In the notation above the 3 index refers to 3-dimensional quantities estimated from the spatial distribution of halos in the UniverseMachine simulation box.The adopted 3 phase-space diagram is independent of projection effects that are inevitable when constructing light cone realisations assuming different observer positions (see Section 3.2).Following commonly used criteria we define infalling halos/galaxies as those or that simultaneously satisfy the following conditions where  ,NFW corresponds to the escape velocity (e.g.Rhee et al. 2017) of a halo assuming an Navarro-Frenk-White profile (NFW, Navarro et al. 1996) and M halo,peak is the maximum historical mass of the halo.Figure 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 Figure 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  ≈ 1 is much higher than that predicted by our baseline semi-empirical model that uses either the Georgakakis et al. ( 2017 Instead, Muñoz Rodríguez et al. (2023) proposed that a log-normal SAR model with mean specific accretion rate log    = −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  ≈ 1.We therefore choose to use the same SAR model in our analysis and apply it only to the infalling galaxies (red circles) of Figure 11.The impact on the AGN radial distribution of the modified SAR for the infalling galaxies is shown in Figure 12.Relative to our baseline model the mean expected number of X-ray selected AGN slightly increases for the radial ring 2 − 2.5  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  500 .This is because the infall population is evenly distributed between  3 = 0.5 − 3  500 as demonstrated by the top panel of Figure 11.In any case, the increase at the ring 2 − 2.5  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  500 among the lightcone realisations with the modified SAR.We find that 20% 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 Figure 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.

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  ≈ 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  500 in massive clusters of galaxies at  ≈ 1 (Koulouridis & Bartalucci 2019) are not physical but instead driven by projection effects.
We select halos at  ≈ 1 in the simulations with masses similar to the clusters of Koulouridis & Bartalucci (2019) and generate mock observations through different sight lines 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: (i) 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  ∼ 0.2.The normalisation of our HODs decreases toward 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.
(ii) 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  ≈ 0.25.The red dots represent recent infall galaxies with dark matter halo masses that have at least 80% of their maximum historical masses (M halo,peak parameter in UniverseMachine catalogue).These are the halos that we consider as infalling in our analysis.The panel at the top shows the (normalised) radial distribution histogram of the different galaxy populations with the same color coding, black refers to the whole population of galaxies, blue to the galaxies in the recent infall region (i.e.those within the blue shaded area) and red for the infall galaxies which dark matter haloes have at least 80% of their maximum historical masses.
(iii) Turning to the projected radial distribution of X-ray selected AGN in the vicinity of massive clusters at  ≈ 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 (Figure 7).
(iv) Our analysis emphasises 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  500 are predicted by the model.Up to 20% of the realisations show amplitudes similar to the observations of Koulouridis & Bartalucci (2019) for massive clusters of galaxies at  ≈ 1 (see Figure 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%.In our model, however, these overdensities in the projected radial distribution are not physical but stochastic and dominated by interlopers (Figure 9).
(v) Fine tuning our model to favour a higher incidence of mock AGN among galaxies in the infall region of masssive halos has little impact to the predicted projected radial distributions (see Figure 12).This further emphasises the significance of sample variance in interpreting projected AGN radial distributions.

Figure 1 .
Figure1.X-ray AGN number counts in the full band (0.5-7.0 keV) for the two selected clusters as function of radial distance in units of  500 .The expected number of sources in the field have been statistically subtracted from each annulus.

Figure 2 .
Figure 2. Graphical demonstration of the semi-empirical model that produces realistic simulated observations of AGN in massive halos.The upper branch shows the construction of the AGN mocks.The starting point (left panel) are dark matter only cosmological simulations.Dark matter halos are populated with galaxies using abundance matching techniques (panel labelled stellar mass).Supermassive black-hole accretion events are painted of top of these galaxies using empirical relations that describe the probability of a galaxy hosting an AGN (panel labelled AGN).See Section 3.1 for details on the mock catalogue construction.The right panel on the upper branch represents a light cone with a field of view 20 arcmin pointing to a massive halo (i.e.cluster of galaxies) in the mock catalogue as described in Section 3.2.The lower branch shows the derivation of the X-ray selection effects.The starting point are X-ray observations (left panel) of cluster of galaxies presented in Section 2. These are processed to obtain the X-ray sensitivity maps (middle panel) explained in Section 2. The anulii in these panels have radii that are multiples of the quantity 0.5 •  500 .These are the anulii used in the observations (see Figure1) to study the radial distribution of AGN in clusters.The different colours represent different rings.The right panel shows the area curve for the different rings which describe the probability of detecting a source at different radial distances from the centre of the cluster.The colour of each line corresponds to the anulii shown in the sensitivity map panel.Both branches merge in the panel labelled "Simulated observation", which represents a simulated field of view that includes the selection effects derived from observations.

Figure 3 .
Figure 3. Examples of two light cones that intersect different sight-lines through the cosmic web.The left panel shows a 3-dimensional projection of the space.The shaded region represents the boundaries of a MDPL2 box.Black points are dark matter haloes with masses > 10 14.25 M ⊙ /ℎ in the simulation.This limit is chosen to help visualisation.Blue and orange points represent the two different light-cones with all the haloes (irrespective of their mass) within their solid angle.The right panels show the corresponding projections of each of the light cones.Each point on the right set of panels represents a dark matter halo.

Figure 4 .
Figure 4. Sketch that illustrates the step-by-step construction of a light-cone for an observer located underneath the centre of the N-body simulation box.The observer's positions is indicated with a black triangle that lies on the plane indicated with the solid horizontal line.In all panels we show a stack of 2 boxes at different snapshot redshifts, coloured differently for illustration purposes.The red shading is for the higher redshift box (relative to the observer) whereas the blue colour correspond to the lower redshift box.The extent of the shaded regions indicates the size of the boxes.The crosses indicate the centres of each box, which in this example are also used as pivot points (see Section 3.2 for details).The dots within each box correspond to dark matter halos with masses M halo > 10 13.5 M ⊙ /h and have the same colour (blue or red) as that of the box they belong to.The black dashed curves represent the iso-redshift surfaces relative to the observer.These define the split-overlap-surfaces used to select halos from the different boxes (see Section 3.2 for details).The construction of the light-cone proceeds from left to right: (i) First, we offset each box in the vertical axis so that the redshift of its pivot point (cross) relative to the observer equals the snapshot redshift of the box; (ii) a set of split-overlap-surface is defined with respect to the grid of boxes (dashed curves); (iii) the set of split-overlap-surfaces is used to remove dublicate halos in the overlap region of the two boxes (middle panel): below the lower dashed curve only blue halos are kept, whereas between the lower and upper dashed curves only red halos remain; (iv) the field-of-view is applied to the box-slices (right panel), keeping only halos that are within a user defined solid angle.
Figure5.A sketch that illustrates the problem of using arbitrary lines of sight when constructing light cones that are forced to intercept a particular halo.The panels show a stack of 3 boxes at different snapshot redshifts.Each of the boxes is coloured differently (blue, green or red) for illustration purposes.The blue and red shadings correspond to the lowest and highest redshifts of the stack.The position of the observer is shown at the bottom of each stack of boxes with the eyeball graphic.The left panel shows an example of a tilted sight line that is forced to contain the position of a halo marked with the star symbol in the green box.This sight line hits the boundaries of the stack of boxes before reaching the expected maximum redshift   (see Section 3.2.1).Such a light cone is clearly incomplete.The right panel visualises a solution to this issue that is based on the construction of two independent light cones.The first one encompasses the region between the observer and the last redshift slice where the light cone is complete.We refer to this component as the foreground light cone.A second independent light cone is then constructed that extends from the previous complete redshift to   .We refer to this component as background light cone.Finally, both light cones are aligned by matching the lines of sight.This is indicated in the far right panel.The arrow shows the rotation that needs to be applied to the background light cone to align it to the foreground one.

Figure 8 .
Figure 8.Each panel plots the radial distributions of mock X-ray AGN (green circles connected with green solid lines) for each of the 100 individual light-cone realisations that point to the massive halo ( 500 ≈ 7.5 × 10 14 M ⊙ /ℎ) with id=7830644447 in UniverseMachine.The semi-empirical model predictions shown in each panel use the specific accretion rate of Georgakakis et al. (2017) to seed galaxies with AGN.The observed projected radial distribution of X-ray selected AGN of the cluster PLCKG266.6-27.3 is shown with the grey/black squares connected with the solid grey/black lines (see Figure1and Section 2).The light-cone realisations that reproduce the observed peak at the distance of 2.5  500 are highlighted by (i) making the observational data points and connecting lines black, (ii) using bold green characters for the light-cone incremental number at the top of the corresponding panel and (iii) change in the background colour from white to grey.

Figure 9 .
Figure 9. Redshift distribution of mock AGN that lie in the radial ring  = 2 − 2.5  500 .Different colours correspond to each of the 7 randomly selected light-cone realisations of Figure 8 (see legend) that reproduce an excess number of projected counts at the radial distance ring  = 2 − 2.5  500 in agreement with the observations presented in Figure 7.

Figure 10 .
Figure 10.The observed projected radial distribution of X-ray selected AGN (black squares and black solid connecting lines) of the clusters PLCKG266.6-27.3(right column of panels) and SPT-CLJ2146-4633 (left column of panels) are compared with the semi-empirical model predictions (coloured lines and shaded regions).The two cluster observations differ in the total Chandra exposure time, with PLCKG266.6-27.3being deeper (about 200 ks) and SPT-CLJ2146-4633 shallower (about 70 ks).In all panels the green lines and shaded regions represent the model that uses the Georgakakis et al. (2017) SARD, while the orange lines and shaded regions are for models that adopt the Aird et al. (2018) SARD for seeding galaxies with AGN.The solid lines are the average of the 100 realisations, while the shaded regions indicate the 1 scatter at fixed radial bin.The model predictions are constructed for two different massive halos in UniverseMachine.The upper row of panels is for the halo with id=7830644447 and mass  500 ≈ 7.5 × 10 14 M ⊙ /ℎ (same as in Figures7, 8).The lower row of panels is for the halo with id=7793510527 and mass  500 ≈ 5 × 10 14 M ⊙ /ℎ (see Figure1).

Figure 11 .
Figure11.The 3-dimensional phase-space diagram used to identify the infalling galaxy population of the cluster with id=7830644447 in UniverseMachine.Black dots correspond to individual galaxies in UniverseMachine.The blue shaded area marks the recent infall region of the parameter space and is defined by the caustic  3 / 500  3 /  = 0.4 (black solid line, e.g.,Kim et al. 2023) and the escape velocity of the equivalent NFW halo profile (dashed black line).The orange shaded area under the caustic  3 / 500  3 /  = 0.4 is often referred to as ancient infall or first infallers region of the phase-space diagram.The red dots represent recent infall galaxies with dark matter halo masses that have at least 80% of their maximum historical masses (M halo,peak parameter in UniverseMachine catalogue).These are the halos that we consider as infalling in our analysis.The panel at the top shows the (normalised) radial distribution histogram of the different galaxy populations with the same color coding, black refers to the whole population of galaxies, blue to the galaxies in the recent infall region (i.e.those within the blue shaded area) and red for the infall galaxies which dark matter haloes have at least 80% of their maximum historical masses.

Figure 12 .
Figure 12.The X-ray AGN radial distribution.The observations (black squares and sold lines) and the semi-empirical model predictions (lines and shaded regions) are plotted at different cluster centric distances normalized to  500 .Black points connected with the solid black line represent the observed radial distribution of the cluster PLCKG266.6-27.3.The solid green line is the mean semi-empirical model prediction in the case of the Georgakakis et al. (2017) specific accretion rate distribution.The magenta solid line correspond to the semi-empirical model in which the modified SARD described in Section 3.1 is applied to the infalling galaxy population identified in Figure11.The light-green shaded and magenta hatched regions within which the semi-empirical model lines are embedded correspond to the 68% confidence intervals of the mean value.They represent the variance between different lines of sight (see the text for further details).Both semi-empirical model predictions are for the massive halo with id=7830644447 in UniverseMachine with virial mass ∼ 8.1 • 10 14 M ⊙ /ℎ.
Aird et al. (2018)our Semi-empirical model predictions for the AGN HOD at different X-ray luminosity thresholds (indicated in the legend) and redshifts (indicated in the title of each panel).The different models of the specific accretion rate distribution used in our semi-empirical approach are indicated with different line styles.Solid lines correspond toAird et al. (2018)whereas dashed lines are for Georgakakis et al. (2017) (see Section 3.1 for details).At the redshift panel  = 0.25 we also compare the predictions of the model with the observational results of Comparat et al. (2023).The black solid line correspond to the best-fit AGN HOD from this study.The 1 uncertainties are shown with by the grey shaded region.semi-empirical model is consistent with the large scale distribution of AGN in the Universe.The AGN HOD, ⟨ (  )|⟩, is defined as the mean number of AGN brighter than the luminosity   in dark matter halos of given mass, .Because of the different halo types (central or satellites) the HOD is usually expressed as a sum of two terms ⟨ (  )|⟩ = ⟨  (  )|⟩ + ⟨  (  )|⟩, ⟨  (  )|⟩ =   •    , (,   ) () , ⟨  (  )|⟩ =    , (  = ,   )   () , Georgakakis et al. (2017)dial distribution of X-ray selected AGN of the cluster PLCKG266.6-27.3(blackpointsconnectedwith the solid black line) is compared with the semi-analytic model predictions (colored lines and shaded regions).The green and orange solid lines show the mean radial distributions of simulated X-ray AGN assuming theGeorgakakis et al. (2017)andAird et al. (