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

Iván Muñoz Rodríguez1,2, Antonis Georgakakis1, Francesco Shankar2, Viola
Allevato3, Silvia Bonoli4,5, Marcella Brusa6,7, Andrea Lapi8 and Akke
Viitanen9
1Institute for Astronomy and Astrophysics, National Observatory of Athens, V.
Paulou & I. Metaxa 11532, Greece
2Department of Physics and Astronomy, University of Southampton, High-field
SO17 1BJ, UK
3INAF - Astronomical Observatory of Capodimonte, Salita Moiariello 16, 80131
Naples, Italy
4Donostia International Physics Center (DIPC), Manuel Lardizabal Ibilbidea, 4,
San Sebastián, Spain
5Ikerbasque, Basque Foundation for Science, E-48013 Bilbao, Spain
6INAF - Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via
Gobetti 93/3, 40129 Bologna, Italy
7Dipartimento di Fisica e Astronomia "Augusto Righi", Università di Bologna,
via Gobetti 93/2, 40129 Bologna, Italy
8SISSA, Via Bonomea 265, 34136 Trieste, Italy
9Department of Physics, Gustaf Hällströmin katu 2, 00014 University of
Helsinki, Finland E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

This paper explores the role of small-scale environment ($<1$ Mpc) in
modulating accretion events onto supermassive black holes by studying the
incidence of Active Galactic Nuclei (AGN) in massive clusters of galaxies. A
flexible, data-driven semi-empirical model is developed based on a minimal set
of parameters and under the zero order assumption that the incidence of AGN in
galaxies is independent of environment. This is used to predict how the
fraction of X-ray selected AGN among galaxies in massive dark matter halos
($\ga 3\times 10^{14}\,M_{\odot}$) evolves with redshift and reveal tensions
with observations. At high redshift, $z\sim 1.2$, the model underpredicts AGN
fractions, particularly at high X-ray luminosities, $L_{X}(\rm 2-10\,keV)\ga
10^{44}\,erg\,s^{-1}$. At low redshift, $z\sim 0.2$, the model estimates
fractions of moderate luminosity AGN ($L_{X}(\rm 2-10\,keV)\ga
10^{43}\,erg\,s^{-1}$) that are a factor of $2-3$ higher than the
observations. These findings reject the zero order assumption on which the
semi-empirical model hinges and point to a strong and redshift-dependent
influence of the small-scale environment on the growth of black holes. Cluster
of galaxies appear to promote AGN activity relative to the model expectation
at $z\sim 1.2$ and suppress it close to the present day. These trends could be
explained by the increasing gas content of galaxies toward higher redshift
combined with an efficient triggering of AGN at earlier times in galaxies that
fall onto clusters.

###### keywords:

galaxies: haloes – galaxies: active – galaxies: clusters: general – galaxies:
nuclei – quasars: supermassive black holes – galaxies: Seyfert – X-rays:
galaxies: clusters

††pubyear: 2022††pagerange: Cosmic evolution of the incidence of Active
Galactic Nuclei in massive clusters: Simulations versus observations–Cosmic
evolution of the incidence of Active Galactic Nuclei in massive clusters:
Simulations versus observations

##  1 Introduction

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

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

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

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

##  2 Observations

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

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

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

##  3 Methodology

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

###  3.1 Galaxy-halo connection: (Sub-)Halo Abundance matching techniques

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

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

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

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

###  3.2 AGN-galaxy connection: specific accretion rate distributions

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

$\lambda_{\rm{sBHAR}}=\frac{\rm{k}_{\rm{bol}}}{1.26\times 10^{38}\times
0.002}\cdot\left(\frac{L_{X}(2-10\,\rm{keV})}{\rm{erg}\,\rm{s}^{-1}}\right)\cdot\left(\frac{M_{\odot}}{M_{\star}}\right).$
(1)

The above equation assumes a Margorrian-type relation between stellar and
black hole mass $\rm{M}_{BH}=0.002\cdot M_{\star}$ (Marconi & Hunt, 2003), an
AGN bolometric correction
$\rm{k}_{\rm{bol}}=L_{\rm{bol}}/L_{X}(2-10\,\rm{keV})=25$ (Elvis et al., 1994)
and the Eddington luminosity of the black hole $1.26\times
10^{38}\,\rm{erg\,s^{-1}}$. The scaling factors in Equation (1) make
$\lambda_{\rm{sBHAR}}$ resemble an Eddington ratio, i.e. the AGN bolometric
luminosity normalised to the Eddington luminosity of the black hole. It is
emphasised that the multiplicative constants in Equation (1) do not affect our
analysis and the assignment of AGN luminosities to UniverseMachine galaxies.

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

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

Figure 1: Specific accretion rate distributions that describe the probability
of a galaxy hosting an AGN with specific accretion rate $\lambda_{\rm sBHAR}$.
The shaded regions correspond to different observational measurements of
P($\lambda_{\rm sBHAR}$). The purple colour shows the Aird et al. (2018)
result, where different line styles indicate different galaxy types: star-
forming galaxies (solid line) or passive (dash-dotted). The green colour shows
the Georgakakis et al. (2017) constraints on the specific accretion are
distribution. All curves correspond to the redshift interval $z=0.5-1$ and
stellar mass interval $M_{\star}=10^{10.5}-10^{11}M_{\odot}$. The extent of
the shaded regions correspond to the 68% confidence interval around the median
(bold curves).

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

$\log\rm{SFR}_{\rm{cut}}=-8.9+0.76\log\left(\frac{M_{\star}}{M_{\odot}}\right)+2.95\left(1+\it{z}\right).$
(2)

At fixed stellar mass and redshift galaxies with star-formation rate above and
below $\rm{SFR}_{\rm{cut}}$ are considered star-forming and passive
respectively.

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

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

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

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

.

###  3.3 Light-cones

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

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

###  3.4 Selection effects

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

####  3.4.1 Cluster sample

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

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

####  3.4.2 Cluster membership

Although in the simulation box the satellite galaxies associated with a given
parent halo are known, we prefer to follow an observational-motivated approach
for defining cluster membership based on the projected and radial distances
relative to the cluster center. Cluster member candidates are selected to lie
within the projected $R_{\rm vir}$ radius of the corresponding halo. This is
similar to the selection of Martini et al. (2009); Martini et al. (2013). We
choose to use $R_{\rm vir}$ instead of $R_{200}$ since they are similar, but
the last is not present in UniverseMachine dataset. We checked that this
assumption does not affect our final results and conclusions.

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

Figure 3: The light-cone of a massive cluster in UniverseMachine at redshift
$z=1.25$ and halo mass $M_{\rm{vir}}=10^{14.51}M_{\odot}/h$. Large circle
marks the virial radius of the cluster, all symbols correspond to
UniverseMachine galaxies within the light-cone. Different symbols and colors
demonstrate selection effects as described in Section 3.4. Purple crosses, red
empty squares and blue diamonds indicate galaxies that are excluded from the
sample because of the selection effects as indicated in the legend. Green
circles and orange triangles indicate the final sample of galaxies and AGNs
with $L_{X}(2-10\,\rm keV)>10^{43}$ erg/s after applying all the selection
effects.

####  3.4.3 Galaxy/AGN sample selection

Table 1: Selection effects adopted for defining the mock cluster galaxy and
AGN samples.

$z_{\rm box}$ | $M_{\rm vir,lim}$ | $m_{\rm lim}$ | band | $M_{\star,\rm lim}$ | LX,lim
---|---|---|---|---|---
(adim.) | ($M_{\odot}/h$) | (mag) | (adim.) | ($M_{\odot}$) | (erg s-1)
(1) | (2) | (3) | (4) | (5) | (6)
1.25 | $3\cdot 10^{14}$ | 18.5 | IRAC1 | 2.2$\cdot$1010 |  $10^{43}$ and $10^{44}$
0.75 | 3.6$\cdot 10^{14}$ | 23.3 | R band | 2.2$\cdot$1010 |  $10^{43}$ and $10^{44}$
0.2 | 5.7$\cdot 10^{14}$ | 19.1 | R band | 2.2$\cdot$1010 |  $10^{43}$ and $10^{44}$

  * (1) Redshift of each cluster corresponding to those of Martini et al. (2009); Martini et al. (2013), (2) Minimum dark matter halo mass (virial) adopted to define the parent cluster sample at different redshifts, (3) minimum magnitude threshold adopted to define galaxy/AGN cluster members at different redshifts, (4) Photometric band used to define the magnitude threshold, (5) Stellar mass limit used to select galaxy/AGN cluster members at different redshifts and (6) X-ray luminosity limit adopted to define the AGN sample at different redshifts.

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

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

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

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

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

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

##  4 Results

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

###  4.1 Impact of selection effects on the incidence of AGN

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

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

###  4.2 X-ray AGN fractions in clusters

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

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

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

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

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

Table 2: Fraction of X-ray AGN relative to galaxies in simulated massive
clusters for different redshifts.

z | LX,lim | f${}_{\rm G+17}\cdot 10^{-2}$ | f${}_{\rm A+18}\cdot 10^{-2}$
---|---|---|---
SM | Mag | SM | Mag
(1) | (2) | (3) | (4)
1.25 | 1043 | 2.48${}^{+0.10}_{-0.13}$ | 2.43${}^{+0.01}_{-0.12}$ | 3.21${}^{+0.02}_{-0.17}$ | 3.08${}^{+0.02}_{-0.16}$
1044 | 0.23 ${}^{+0.04}_{-0.04}$ | 0.22${}^{+0.04}_{-0.04}$ | 0.28${}^{+0.04}_{-0.04}$ | 0.27 ${}^{+0.04}_{-0.04}$
0.75 | 1043 | 1.56 ${}^{+0.03}_{-0.03}$ | 1.04${}^{+0.020}_{-0.021}$ | 1.21${}^{+0.03}_{-0.03}$ | 0.635${}^{+0.020}_{-0.016}$
1044 | 0.096${}^{+0.007}_{-0.007}$ | 0.063${}^{+0.005}_{-0.005}$ | 0.062${}^{+0.006}_{-0.006}$ | 0.038${}^{+0.004}_{-0.004}$
0.2 | 1043 | 0.300${}^{+0.007}_{-0.007}$ | 0.266${}^{+0.006}_{-0.006}$ | 0.328${}^{+0.008}_{-0.008}$ | 0.208${}^{+0.006}_{-0.006}$
1044 | 0.084${}^{+0.004}_{-0.004}$ | 0.072${}^{+0.003}_{-0.003}$ | 0.036${}^{+0.003}_{-0.003}$ | 0.0263${}^{+0.0020}_{-0.0020}$

  * (1) Redshift, (2) X-ray luminosity threshold on the $2-10$ keV band, adopted for the AGN sample, in units of erg/s, (3) fractions corresponding to the Georgakakis et al. (2017) model, for the stellar mass (SM) and magnitude (Mag) selected samples, (4) fractions corresponding to Aird et al. (2018) model, for the stellar mass (SM) and magnitude (Mag) selected samples.

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

##  5 Discussion

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

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

###  5.1 AGN fractions in high redshift clusters

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

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

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

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

###  5.2 Reproducing high-$z$ AGN fractions in clusters

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

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

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

###  5.3 Evolution of AGN incidence in cluster vs. field

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

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

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

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

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

###  5.4 Physical interpretation

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

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

  1. 1.

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

  2. 2.

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

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

##  6 Summary and conclusions

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

  1. 1.

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

  2. 2.

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

  3. 3.

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

  4. 4.

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

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

## Acknowledgements

The authors thank the anonymous referee for their careful reading and their
comments of the paper. We are grateful as well to Ángel Ruiz Camuñas for his
helpful comments and advises. 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. AL is partly supported by the PRIN
MIUR 2017 prot. 20173ML3WW 002 ‘Opening the ALMA window on the cosmic
evolution of gas, stars, and massive black holes’. SB acknowledges the project
PGC2018-097585-B-C22, MINECO/FEDER, UE of the Spanish Ministerio de Economia,
Industria y Competitividad. This research made use of
Astropy,444http://www.astropy.org a community-developed core Python package
for Astronomy (Astropy Collaboration et al., 2013, 2018);
Colossus555https://bdiemer.bitbucket.io/colossus/index.html (Diemer, 2018),
Numpy666https://numpy.org/ (van der Walt et al., 2011),
Scipy777https://scipy.org/ (Virtanen et al., 2020),
Matplotlib888https://matplotlib.org/ (Hunter, 2007) and
HaloMod999https://pypi.org/project/halomod/ (Murray et al., 2013, 2021). For
the purpose of open access, the authors have applied a CC-BY public copyright
license to any author accepted manuscript version arising.

## Data Availability

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

## References

  * Aird & Coil (2021) Aird J., Coil A. L., 2021, MNRAS, 502, 5962
  * Aird et al. (2012) Aird J., et al., 2012, ApJ, 746, 90
  * Aird et al. (2015) Aird J., Coil A. L., Georgakakis A., Nandra K., Barro G., Pérez-González P. G., 2015, MNRAS, 451, 1892
  * Aird et al. (2018) Aird J., Coil A. L., Georgakakis A., 2018, MNRAS, 474, 1225
  * Alberts et al. (2016) Alberts S., et al., 2016, ApJ, 825, 72
  * Alberts et al. (2022) Alberts S., Adams J., Gregg B., Pope A., Williams C. C., Eisenhardt P. R. M., 2022, ApJ, 927, 235
  * Allevato et al. (2012) Allevato V., et al., 2012, ApJ, 758, 47
  * Allevato et al. (2021) Allevato V., Shankar F., Marsden C., Rasulov U., Viitanen A., Georgakakis A., Ferrara A., Finoguenov A., 2021, ApJ, 916, 34
  * Arthur et al. (2019) Arthur J., et al., 2019, MNRAS, 484, 3968
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Becker et al. (2007) Becker M. R., et al., 2007, ApJ, 669, 905
  * Behroozi et al. (2010) Behroozi P. S., Conroy C., Wechsler R. H., 2010, ApJ, 717, 379
  * Behroozi et al. (2013a) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013a, ApJ, 762, 109
  * Behroozi et al. (2013b) Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A., Primack J. R., 2013b, ApJ, 763, 18
  * Behroozi et al. (2015) Behroozi P. S., et al., 2015, Monthly Notices of the Royal Astronomical Society, 450, 1546
  * Behroozi et al. (2019) Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143
  * Bongiorno et al. (2012) Bongiorno A., et al., 2012, MNRAS, 427, 3103
  * Bournaud et al. (2012) Bournaud F., et al., 2012, ApJ, 757, 81
  * Brandt & Alexander (2015) Brandt W. N., Alexander D. M., 2015, A&ARv, 23, 1
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Bufanda et al. (2017) Bufanda E., et al., 2017, MNRAS, 465, 2531
  * Butcher & Oemler (1984) Butcher H., Oemler A. J., 1984, ApJ, 285, 426
  * Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, ApJ, 533, 682
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Christlein & Zabludoff (2003) Christlein D., Zabludoff A. I., 2003, ApJ, 591, 764
  * Chung et al. (2011) Chung S. M., Eisenhardt P. R., Gonzalez A. H., Stanford S. A., Brodwin M., Stern D., Jarrett T., 2011, ApJ, 743, 34
  * Ciotti & Ostriker (2007) Ciotti L., Ostriker J. P., 2007, ApJ, 665, 1038
  * Cisternas et al. (2015) Cisternas M., Sheth K., Salvato M., Knapen J. H., Civano F., Santini P., 2015, ApJ, 802, 137
  * DeGraf et al. (2017) DeGraf C., Dekel A., Gabor J., Bournaud F., 2017, MNRAS, 466, 1462
  * Dekel & Birnboim (2006) Dekel A., Birnboim Y., 2006, MNRAS, 368, 2
  * Di Matteo et al. (2005) Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
  * Diemer (2018) Diemer B., 2018, ApJS, 239, 35
  * Digby-North et al. (2010) Digby-North J. A., et al., 2010, MNRAS, 407, 846
  * Eastman et al. (2007) Eastman J., Martini P., Sivakoff G., Kelson D. D., Mulchaey J. S., Tran K.-V., 2007, ApJ, 664, L9
  * Eisenhardt et al. (2008) Eisenhardt P. R. M., et al., 2008, ApJ, 684, 905
  * Elbaz et al. (2007) Elbaz D., et al., 2007, A&A, 468, 33
  * Elvis et al. (1994) Elvis M., et al., 1994, ApJS, 95, 1
  * Fabian (1994) Fabian A. C., 1994, ARA&A, 32, 277
  * Fakhouri & Ma (2008) Fakhouri O., Ma C.-P., 2008, MNRAS, 386, 577
  * Gatti et al. (2016) Gatti M., Shankar F., Bouillot V., Menci N., Lamastra A., Hirschmann M., Fiore F., 2016, MNRAS, 456, 1073
  * Georgakakis et al. (2011) Georgakakis A., et al., 2011, MNRAS, 418, 2590
  * Georgakakis et al. (2017) Georgakakis A., Aird J., Schulze A., Dwelly T., Salvato M., Nandra K., Merloni A., Schneider D. P., 2017, MNRAS, 471, 1976
  * Georgakakis et al. (2019) Georgakakis A., Comparat J., Merloni A., Ciesla L., Aird J., Finoguenov A., 2019, MNRAS, 487, 275
  * Georgakakis et al. (2020) Georgakakis A., Ruiz A., LaMassa S. M., 2020, MNRAS, 499, 710
  * Gobat et al. (2020) Gobat R., Magdis G., D’Eugenio C., Valentino F., 2020, A&A, 644, L7
  * Gottlöber et al. (2001) Gottlöber S., Klypin A., Kravtsov A. V., 2001, ApJ, 546, 223
  * Goulding et al. (2014) Goulding A. D., et al., 2014, ApJ, 783, 40
  * Haggard et al. (2010) Haggard D., Green P. J., Anderson S. F., Constantin A., Aldcroft T. L., Kim D.-W., Barkhouse W. A., 2010, The Astrophysical Journal, 723, 1447
  * Haines et al. (2012) Haines C. P., et al., 2012, ApJ, 754, 97
  * Hopkins & Hernquist (2006) Hopkins P. F., Hernquist L., 2006, ApJS, 166, 1
  * Hunter (2007) Hunter J. D., 2007, Computing in Science Engineering, 9, 90
  * Ilbert et al. (2009) Ilbert O., et al., 2009, ApJ, 690, 1236
  * Ilbert et al. (2013) Ilbert O., et al., 2013, A&A, 556, A55
  * Kauffmann & Heckman (2009) Kauffmann G., Heckman T. M., 2009, MNRAS, 397, 135
  * Kauffmann et al. (2004) Kauffmann G., White S. D. M., Heckman T. M., Ménard B., Brinchmann J., Charlot S., Tremonti C., Brinkmann J., 2004, MNRAS, 353, 713
  * Kereš et al. (2005) Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2
  * Kereš et al. (2009) Kereš D., Katz N., Fardal M., Davé R., Weinberg D. H., 2009, MNRAS, 395, 160
  * Klypin et al. (2011) Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102
  * Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Kotecha et al. (2022) Kotecha S., et al., 2022, MNRAS, 512, 926
  * Koulouridis & Bartalucci (2019) Koulouridis E., Bartalucci I., 2019, A&A, 623, L10
  * Koulouridis & Plionis (2010) Koulouridis E., Plionis M., 2010, ApJ, 714, L181
  * Kravtsov et al. (2004) Kravtsov A. V., Berlind A. A., Wechsler R. H., Klypin A. A., Gottlöber S., Allgood B., Primack J. R., 2004, ApJ, 609, 35
  * Krishnan et al. (2017) Krishnan C., et al., 2017, MNRAS, 470, 2170
  * Lapi et al. (2011) Lapi A., et al., 2011, ApJ, 742, 24
  * Lapi et al. (2018) Lapi A., et al., 2018, ApJ, 857, 22
  * Lehmer et al. (2009) Lehmer B. D., et al., 2009, ApJ, 691, 687
  * Lemson & Virgo Consortium (2006) Lemson G., Virgo Consortium t., 2006, arXiv e-prints, pp astro–ph/0608019
  * Lopes et al. (2017) Lopes P. A. A., Ribeiro A. L. B., Rembold S. B., 2017, MNRAS, 472, 409
  * Man et al. (2019) Man Z.-y., Peng Y.-j., Kong X., Guo K.-x., Zhang C.-p., Dou J., 2019, MNRAS, 488, 89
  * Mancone et al. (2010) Mancone C. L., Gonzalez A. H., Brodwin M., Stanford S. A., Eisenhardt P. R. M., Stern D., Jones C., 2010, ApJ, 720, 284
  * Marconi & Hunt (2003) Marconi A., Hunt L. K., 2003, ApJ, 589, L21
  * Marconi et al. (2004) Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, MNRAS, 351, 169
  * Marshall et al. (2018) Marshall M. A., Shabala S. S., Krause M. G. H., Pimbblet K. A., Croton D. J., Owers M. S., 2018, MNRAS, 474, 3615
  * Martini et al. (2006) Martini P., Kelson D. D., Kim E., Mulchaey J. S., Athey A. A., 2006, ApJ, 644, 116
  * Martini et al. (2009) Martini P., Sivakoff G. R., Mulchaey J. S., 2009, ApJ, 701, 66
  * Martini et al. (2013) Martini P., et al., 2013, ApJ, 768, 1
  * McBride et al. (2009) McBride J., Fakhouri O., Ma C.-P., 2009, MNRAS, 398, 1858
  * Merloni & Heinz (2008) Merloni A., Heinz S., 2008, MNRAS, 388, 1011
  * Miller et al. (2003) Miller C. J., Nichol R. C., Gómez P. L., Hopkins A. M., Bernardi M., 2003, ApJ, 597, 142
  * Mishra & Dai (2020) Mishra H. D., Dai X., 2020, AJ, 159, 69
  * Miyaji et al. (2011) Miyaji T., Krumpe M., Coil A. L., Aceves H., 2011, ApJ, 726, 83
  * Moster et al. (2010) Moster B. P., Somerville R. S., Maulbetsch C., van den Bosch F. C., Macciò A. V., Naab T., Oser L., 2010, ApJ, 710, 903
  * Moster et al. (2018) Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822
  * Moustakas et al. (2013) Moustakas J., et al., 2013, The Astrophysical Journal, 767, 50
  * Munari et al. (2013) Munari E., Biviano A., Borgani S., Murante G., Fabjan D., 2013, MNRAS, 430, 2638
  * Murray et al. (2013) Murray S. G., Power C., Robotham A. S. G., 2013, Astronomy and Computing, 3, 23
  * Murray et al. (2021) Murray S. G., Diemer B., Chen Z., Neuhold A. G., Schnapp M. A., Peruzzi T., Blevins D., Engelman T., 2021, Astronomy and Computing, 36, 100487
  * Muzzin et al. (2013) Muzzin A., et al., 2013, ApJ, 777, 18
  * Nantais et al. (2016) Nantais J. B., et al., 2016, A&A, 592, A161
  * Noble et al. (2019) Noble A. G., et al., 2019, ApJ, 870, 56
  * Padovani et al. (2017) Padovani P., et al., 2017, A&ARv, 25, 2
  * Papovich et al. (2018) Papovich C., et al., 2018, ApJ, 854, 30
  * Pearl et al. (2021) Pearl A. N., Bezanson R., Zentner A. R., Newman J. A., Goulding A. D., Whitaker K. E., Johnson S. D., Greene J. E., 2021, arXiv e-prints, p. arXiv:2112.00035
  * Peluso et al. (2022) Peluso G., et al., 2022, ApJ, 927, 130
  * Pimbblet et al. (2013) Pimbblet K. A., Shabala S. S., Haines C. P., Fraser-McKelvie A., Floyd D. J. E., 2013, MNRAS, 429, 1827
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Poggianti et al. (2017) Poggianti B. M., et al., 2017, Nature, 548, 304
  * Popesso & Biviano (2006) Popesso P., Biviano A., 2006, A&A, 460, L23
  * Popesso et al. (2015) Popesso P., et al., 2015, A&A, 579, A132
  * Ricarte et al. (2020) Ricarte A., Tremmel M., Natarajan P., Quinn T., 2020, ApJ, 895, L8
  * Richardson et al. (2012) Richardson J., Zheng Z., Chatterjee S., Nagai D., Shen Y., 2012, ApJ, 755, 30
  * Saintonge et al. (2013) Saintonge A., et al., 2013, ApJ, 778, 2
  * Santini et al. (2014) Santini P., et al., 2014, A&A, 562, A30
  * Schreiber et al. (2015) Schreiber C., et al., 2015, A&A, 575, A74
  * Schulze et al. (2015) Schulze A., et al., 2015, MNRAS, 447, 2085
  * Shankar et al. (2020) Shankar F., et al., 2020, Nature Astronomy, 4, 282
  * Shen et al. (2013) Shen Y., et al., 2013, ApJ, 778, 98
  * Soltan (1982) Soltan A., 1982, MNRAS, 200, 115
  * Tacconi et al. (2010) Tacconi L. J., et al., 2010, Nature, 463, 781
  * Tasitsiomi et al. (2004) Tasitsiomi A., Kravtsov A. V., Wechsler R. H., Primack J. R., 2004, ApJ, 614, 533
  * Ueda et al. (2014) Ueda Y., Akiyama M., Hasinger G., Miyaji T., Watson M. G., 2014, ApJ, 786, 104
  * Vale & Ostriker (2004) Vale A., Ostriker J. P., 2004, MNRAS, 353, 189
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vulcani et al. (2013) Vulcani B., et al., 2013, A&A, 550, A58
  * Vulcani et al. (2018) Vulcani B., et al., 2018, MNRAS, 480, 3152
  * Webb et al. (2013) Webb T. M. A., et al., 2013, AJ, 146, 84
  * Wetzel et al. (2013) Wetzel A. R., Tinker J. L., Conroy C., van den Bosch F. C., 2013, MNRAS, 432, 336
  * Williams et al. (2022) Williams C. C., et al., 2022, ApJ, 929, 35
  * Yang et al. (2018) Yang L., Tozzi P., Yu H., Lusso E., Gaspari M., Gilli R., Nardini E., Risaliti G., 2018, ApJ, 859, 65
  * Zinger et al. (2018) Zinger E., Dekel A., Kravtsov A. V., Nagai D., 2018, MNRAS, 475, 3654
  * de Souza et al. (2016) de Souza R. S., et al., 2016, MNRAS, 461, 2115
  * van der Burg et al. (2020) van der Burg R. F. J., et al., 2020, A&A, 638, A112
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science Engineering, 13, 22
  * von der Linden et al. (2010) von der Linden A., Wild V., Kauffmann G., White S. D. M., Weinmann S., 2010, MNRAS, 404, 1231