Draft version April 2, 2024

Typeset using LATEX twocolumn style in AASTeX63

Mapping the Growth of Supermassive Black Holes as a Function of Galaxy Stellar Mass and Redshift
Fan Zou,1, 2 Zhibo Yu,1, 2 W. N. Brandt,1, 2, 3 Hyungsuk Tak,4, 1, 5 Guang Yang,6, 7 and Qingling Ni8
1 Department of Astronomy and Astrophysics, 525 Davey Lab, The Pennsylvania State University, University Park, PA 16802, USA
2 Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA
3 Department of Physics, 104 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA
arXiv:2404.00097v1 [astro-ph.GA] 29 Mar 2024

4 Department of Statistics, The Pennsylvania State University, University Park, PA 16802, USA
5 Institute for Computational and Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA
6 Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands
7 SRON Netherlands Institute for Space Research, Postbus 800, 9700 AV Groningen, The Netherlands
8 Max-Planck-Institut für extraterrestrische Physik (MPE), Gießenbachstraße 1, D-85748 Garching bei München, Germany

ABSTRACT
The growth of supermassive black holes is strongly linked to their galaxies. It has been shown that the
population mean black-hole accretion rate (BHAR) primarily correlates with the galaxy stellar mass (M⋆ ) and
redshift for the general galaxy population. This work aims to provide the best measurements of BHAR as a
function of M⋆ and redshift over ranges of 109.5 < M⋆ < 1012 M⊙ and z < 4. We compile an unprecedentedly
large sample with eight thousand active galactic nuclei (AGNs) and 1.3 million normal galaxies from nine high-
quality survey fields following a wedding-cake design. We further develop a semiparametric Bayesian method
that can reasonably estimate BHAR and the corresponding uncertainties, even for sparsely populated regions in
the parameter space. BHAR is constrained by X-ray surveys sampling the AGN accretion power and UV-to-
infrared multi-wavelength surveys sampling the galaxy population. Our results can independently predict the
X-ray luminosity function (XLF) from the galaxy stellar mass function (SMF), and the prediction is consistent
with the observed XLF. We also try adding external constraints from the observed SMF and XLF. We further
measure BHAR for star-forming and quiescent galaxies and show that star-forming BHAR is generally larger
than or at least comparable to the quiescent BHAR.

Keywords: Supermassive black holes (1663); X-ray active galactic nuclei (2035); Galaxies (573)

1. INTRODUCTION Alexander 2015; Brandt & Yang 2022). Therefore, X-ray

Supermassive black holes (SMBHs) and galaxies appear to surveys can be used to constrain the accretion distribution
be fundamentally linked (e.g., Kormendy & Ho 2013). Es- and the black-hole accretion rate (BHAR = dMBH /dt) of
pecially, the SMBH mass (MBH ) and the galaxy bulge mass SMBHs (e.g., Aird et al. 2012, 2018; Bongiorno et al. 2012,
are found to be tightly correlated in the local universe, and 2016; Georgakakis et al. 2017; Wang et al. 2017; Yang et al.
the cosmic evolution of the SMBH accretion density and the 2017, 2018, 2019; Ni et al. 2019, 2021a). However, the dura-
star-formation density are broadly similar and both peak at tion of a galaxy within the AGN phase is relatively short, and
z ≈ 2. Therefore, it is a fundamental question in extragalactic AGNs also have strong variability (e.g., Hickox et al. 2014;
astronomy to understand how cosmic SMBH growth depends Yang et al. 2016); thus, BHAR is often averaged over a large
on galaxy properties. galaxy sample as an estimator of BHAR, the long-term pop-
SMBHs grow primarily through rapid accretion when they ulation mean BHAR.
are observed as active galactic nuclei (AGNs); mergers are BHAR has been shown to be redshift-dependent and cor-
an additional growth mode. X-ray emission is known to be related with both stellar mass (M⋆ ) and star-formation rate
a good indicator of AGN activity because of its universal- (SFR), while the M⋆ dependence is more fundamental for
ity among AGNs, high penetrating power through obscura- the general galaxy population (e.g., Hickox et al. 2014; Yang
tion, and low dilution from galaxy starlight (e.g., Brandt & et al. 2017, 2018). In recent years, it has been found that
BHAR is also strongly related to galaxy morphology (e.g.,
Ni et al. 2019, 2021a; Yang et al. 2019; Aird et al. 2022),
E-mail: [email protected] which may be more fundamental than the BHAR − M⋆ re-
2

lation. However, such morphological measurements are ex- 2. DATA AND SAMPLE
pensive to obtain and require superb image resolution from, We use the data within the Cosmic Assembly Near-
e.g., the Hubble Space Telescope, which inevitably are lim- Infrared Deep Extragalactic Legacy Survey (CANDELS)
ited to small sky areas and thus can only provide a limited fields, four of the Vera C. Rubin Observatory Legacy Sur-
sample size covering a limited parameter space. In con- vey of Space and Time (LSST) Deep-Drilling Fields (DDFs),
trast, M⋆ and SFR are much more accessible. Notably, mod- and eROSITA Final Equatorial Depth Survey (eFEDS) field.
ern multi-wavelength photometric surveys have provided ex- CANDELS and the LSST DDFs both contain several distinct
tensive photometric data, based on which M⋆ and SFR can fields, and we put those individual fields sharing similar areas
be estimated by fitting the corresponding spectral energy and depths under the same umbrella (“CANDELS” or “LSST
distributions (SEDs; e.g., Zou et al. 2022). Therefore, a DDFs”) for convenience. Our adopted fields have X-ray
well-measured BHAR − M⋆ relation is still essential to link coverage to provide AGN information and quality multi-
SMBHs and galaxies. wavelength data, which are essential for measuring galaxy
The latest version of BHAR as a function of (M⋆ , z) properties. We summarize our field information in Table 1
was derived in Yang et al. (2018) using the data from the and discuss them in the following subsections.
Cosmic Evolution Survey (COSMOS; Laigle et al. 2016;
Weaver et al. 2022), Great Observatories Origins Deep Sur- 2.1. CANDELS Fields
vey (GOODS)-S, and GOODS-N (Grogin et al. 2011; Koeke-
CANDELS (Grogin et al. 2011; Koekemoer et al. 2011)
moer et al. 2011). Although the relation in Yang et al. (2018)
is a pencil-beam survey covering five fields – GOODS-S
has proved to be successful over the years, there are still
(0.05 deg2 ), GOODS-N (0.05 deg2 ), Extended Groth Strip
some limitations. First, although COSMOS, GOODS-S, and
(EGS; 0.06 deg2 ), UKIRT Infrared Deep Sky Survey Ultra-
GOODS-N are sufficiently deep to probe BHAR up to z ≈ 4,
Deep Survey (UDS; 0.06 deg2 ), and a tiny part of COSMOS
they cannot effectively sample the last half of cosmic time
(denoted as CANDELS-COSMOS hereafter; 0.06 deg2 ). All
(z ≲ 0.8) because of their limited sky areas. Second, Yang
the fields have ultra-deep UV-to-IR data (see, e.g., Yang
et al. (2018) adopted strong assumptions when parametri-
et al. 2019 and references therein), allowing for detections
cally estimating BHAR, which may lead to underestimated
of galaxies up to high redshift and low M⋆ and reliable
BHAR uncertainties. Built upon Yang et al. (2018), this work
measurements of these galaxies’ properties. The first four
aims to provide the best available map of BHAR(M⋆ , z) with
have deep Chandra observations reaching ∼Ms depths from
the currently best data and statistical methodology. Now is
Luo et al. (2017; GOODS-S), Xue et al. (2016; GOODS-
indeed the right time to re-measure BHAR given the fact
N), Nandra et al. (2015; EGS), and Kocevski et al. (2018;
that the past five years have witnessed the completion of
UDS) and can thus effectively sample AGNs at high red-
several large, sensitive X-ray surveys in fields together with
shift and/or low luminosity. However, CANDELS-COSMOS
deep optical-to-IR surveys (e.g., Chen et al. 2018; Ni et al.
shares the same X-ray depth as the full COSMOS field, and
2021b). These new X-ray surveys, when combined with pre-
CANDELS-COSMOS is much smaller. Therefore, we do not
vious ones, can return a large AGN sample over ten times
use CANDELS-COSMOS but will instead directly analyze
larger than previous ones, as will be discussed in Section 2.
the full COSMOS field in Section 2.2.
In this work, we compile an unprecedentedly large sample
We adopt the galaxy-property catalog in Yang et al. (2019),
from nine well-studied survey fields, many of which were
who have measured M⋆ and SFRs by fitting SEDs for all the
finished after Yang et al. (2018) and even within ≲ 2 yrs
CANDELS sources.
before this work. Our adopted surveys follow a wedding-
cake design and contain both deep, pencil-beam and shallow, 2.2. LSST DDFs
wide ones, allowing us to effectively explore a wide range
The LSST DDFs (e.g., Brandt et al. 2018; Zou et al. 2022)
of parameter space. We further develop a semiparametric
include five fields – COSMOS, Wide Chandra Deep Field-
Bayesian approach that can return reasonable estimations and
South (W-CDF-S), European Large-Area Infrared Space
uncertainties, even for sparsely populated regions in the pa-
Observatory Survey-S1 (ELAIS-S1), XMM-Newton Large
rameter space.
Scale Structure (XMM-LSS), and Euclid Deep Field-South
This work is structured as follows. Section 2 describes
(EDF-S). EDF-S has been selected as a LSST DDF only re-
the data. Section 3 presents our methodology and BHAR
cently in 2022 and currently does not have sufficient data
measurements. In Section 4, we discuss the implications of
available, and we thus only use the former four original LSST
our results. Section 5 summarizes this work. We adopt a flat
DDFs with superb data accumulated over ∼ a decade. Note
ΛCDM cosmology with H0 = 70 km s−1 Mpc−1 , ΩΛ = 0.7,
that this work does not use any actual LSST data because the
and Ω M = 0.3.
Vera C. Rubin Observatory is still under construction at the
time of writing this article.
 Table 1. Basic Information for the Fields Used in This Work
Field Area mlim X-ray Depth X-ray ref. Galaxy ref. Photo-z ref. AGN Galaxies (a, b)
(deg2 ) (AB mag) (ks)
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
GOODS-S 0.05 26.5 (H) 7000 (Chandra) 5 1 4,8 224 (111) 4144 (−15.87, 2.63)
GOODS-N 0.05 26.5 (H) 2000 (Chandra) 8 1 1,11 174 (167) 4603 (−15.49, 2.58)
EGS 0.06 26.5 (H) 800 (Chandra) 6 1 9 112 (10) 5889 (−15.13, 3.08)
UDS 0.06 26.5 (H) 600 (Chandra) 4 1 8 117 (25) 5010 (−15.05, 4.90)
COSMOS 1.27 24 (K s ) 160 (Chandra) 3 2 5,10 1459 (880) 86765 (−14.68, 5.19)
ELAIS-S1 2.93 23.5 (K s ) 30 (XMM-Newton) 7 3 6,12 676 (261) 157791 (−13.90, 4.57)
W-CDF-S 4.23 23.5 (K s ) 30 (XMM-Newton) 7 3 6,12 872 (311) 210727 (−13.86, 4.97)
XMM-LSS 4.20 23.5 (K s ) 40 (XMM-Newton) 2 3 2 1765 (898) 254687 (−14.09, 5.36)
eFEDS 59.75 22 (Z) 2 (eROSITA) 1 2 3,7 2667 (1156) 615068 (−13.51, 2.59)
Notes. (1) Field names. GOODS-S, GOODS-N, EGS, and UDS belong to CANDELS and are discussed in Section 2.1. COSMOS, ELAIS-S1, W-CDF-S,
and XMM-LSS belong to the LSST DDFs and are discussed in Section 2.2. eFEDS is discussed in Section 2.3. (2) Sky areas of the fields, only accounting
for the regions we are using. (3) The limiting AB magnitudes we adopted in Section 2.4 to calculate the M⋆ completeness curves, and the reference bands
are written within parentheses. (4) The typical depths in exposure time of the X-ray surveys, and the parentheses list the observatories with which our
adopted X-ray surveys were conducted. For XMM-Newton, the reported exposure is the typical flare-filtered one for a single EPIC camera. All the three
EPIC cameras (one EPIC-pn and two EPIC-MOS) were used for the XMM-Newton observations, adding to ≈ 80 − 100 ks EPIC exposure in total. The
“a” parameter values in Column 10 of this table represent typical flux limits in 2 − 10 keV. (5) The references for the X-ray surveys. (6) The references
for our adopted host-galaxy properties. All of these references have appropriately considered the AGN emission for AGNs. (7) Representative references
examining the photo-zs in the corresponding fields. (8) Numbers of AGNs. The parentheses list the numbers of sources with spec-zs. The surface number
density of eFEDS AGNs is much smaller than those in the other fields primarily because the eFEDS limiting magnitude is much brighter. (9) Numbers of
normal galaxies. (10) The parameters describing the X-ray detection function; see Equation 3. There is a subtle difference between eFEDS and other fields
– the eFEDS X-ray detection function is for the intrinsic 2 − 10 keV flux, while the others’ are for the observed flux.
X-ray references. (1) Brunner et al. (2022); (2) Chen et al. (2018); (3) Civano et al. (2016); (4) Kocevski et al. (2018); (5) Luo et al. (2017); (6) Nandra
et al. (2015); (7) Ni et al. (2021b); (8) Xue et al. (2016).
Galaxy references. (1) Yang et al. (2019); (2) Yu et al. (2023); (3) Zou et al. (2022).
Photo-z references. (1) Barro et al. (2019); (2) Chen et al. (2018); (3) Driver et al. (2022); (4) Luo et al. (2017); (5) Marchesi et al. (2016); (6) Ni et al.
(2021b); (7) Salvato et al. (2022); (8) Santini et al. (2015); (9) Stefanon et al. (2017); (10) Weaver et al. (2022); (11) Xue et al. (2016); (12) Zou et al.
(2021).
3
4

COSMOS is a deg2 -scale field with deep multi-wavelength galaxies (AGNs).1 Spec-zs are adopted when available, and
data (e.g., Weaver et al. 2022). Civano et al. (2016) pre- half of the involved AGNs have spec-zs.
sented medium-depth (≈ 160 ks) Chandra observations in We select sources with 0.05 < z < 4 and
log M⋆,min = 9.5 < log M⋆ < log M⋆,max = 12 . Sources



COSMOS. The galaxy properties measured through SED fit-
ting covering the X-ray to far-IR are taken from Yu et al. labeled as stars are removed, as has been presented in the
(2023). We only use the region with “FLAG COMBINED references in Column 6 of Table 1. Only ≲ 15% of sources
= 0” (i.e., within the UltraVISTA region and far from bright in each field are classified as stars. We apply a lower cut for
stars and image edges) in Weaver et al. (2022) to ensure qual- z because photo-zs are less reliable when too small (e.g., see
ity multi-wavelength characterizations. Ni et al. (2021b) pre- Appendix C of Zou et al. 2021), and the peculiar motions be-
sented ≈ 30 ks XMM-Newton observations in ELAIS-S1 and come non-negligible as well. We limit log M⋆ > 9.5 because
W-CDF-S, and Chen et al. (2018) presented ≈ 40 ks XMM- dwarf AGNs usually have much less reliable measurements
Newton observations in XMM-LSS. The galaxy properties in and require special treatment (e.g., Zou et al. 2023). We ap-
these three fields are taken from Zou et al. (2022). We limit ply the same upper cuts as in Yang et al. (2018) for both M⋆
our analyses to the overlapping region between the X-ray cat- and z because very few sources can exceed these thresholds.
alogs and Zou et al. (2022) because quality multi-wavelength We further construct a complete sample by applying
data are essential for estimating photometric redshifts (photo- redshift-dependent M⋆ cuts. To estimate the M⋆ depth for
zs), M⋆ , and SFRs. Besides, GOODS-S and UDS in Sec- each field, we first adopt a reference band and denote its lim-
tion 2.1 are embedded within W-CDF-S and XMM-LSS, re- iting magnitude as mlim . Following Pozzetti et al. (2010),
spectively, and we remove the regions covered by GOODS-S we convert the magnitude depth to the expected limiting
and UDS to avoid double counting sources. Due to these M⋆ for each galaxy with a magnitude of m: log Mlim =
reasons, our areas are slightly smaller than those reported in log M⋆ + 0.4(m − mlim ). At each redshift, we adopt the M⋆
Chen et al. (2018) and Ni et al. (2021b). completeness threshold as the value above which 90% of the
Mlim values lie. Sources below the M⋆ completeness curves
are removed. For the CANDELS fields, we adopt the H band
with a limiting magnitude of 26.5 mag, and almost all the
2.3. eFEDS sources above our log M⋆ cut of 9.5 are above the CANDELS
2 2
eFEDS is a 10 deg -scale field covered by eROSITA with M⋆ completeness curves, enabling constraints upon BHAR
≈ 2 ks observations (Brunner et al. 2022). We focus on in the low-M⋆ and high-z regime. For the LSST DDFs, we
the 60 deg2 GAMA09 region (Driver et al. 2022) within adopt the K s band, and their limiting K s magnitudes are 24
eFEDS because the remaining parts do not have sufficient for COSMOS (Laigle et al. 2016) and 23.5 for W-CDF-S,
multi-wavelength data to constrain the host-galaxy properties ELAIS-S1, and XMM-LSS (Jarvis et al. 2013), respectively.
(e.g., Salvato et al. 2022). Unlike Chandra or XMM-Newton, For eFEDS, we adopt the Z band with a limiting magnitude
eROSITA mostly works at < 2.3 keV, which is more sen- of 22. These M⋆ completeness cuts also automatically ensure
sitive to obscuration. We thus rely on the X-ray properties the general SED-fitting reliability. The typical i-band magni-
cataloged in Liu et al. (2022) for eFEDS sources, which are tudes of sources at these limiting magnitudes are i ≈ 24.8 at
measured through detailed X-ray spectral fitting and thereby K s = 23.5, i ≈ 25.3 at K s = 24, and i ≈ 22.4 at Z = 22.
can largely overcome obscuration effects. As suggested in These i-band magnitudes are roughly equal to the nominal
Liu et al. (2022), we only use sources with detection likeli- “depths” of SEDs in Zou et al. (2022; see their Figure 30)
hoods > 10 because fainter sources generally do not allow and Yu et al. (2023), below which the number of available
meaningful X-ray spectral analyses. photometric bands may become small.
We then define λ = LX /M⋆ , where LX is the intrinsic
2 − 10 keV luminosity, and we always adopt erg s−1 M⊙−1 as
the unit for λ. We use the X-ray surveys mentioned in the
2.4. Sample Construction previous subsections to select AGNs. Following Aird et al.
(2012) and Yang et al. (2018), we only use sources detected
Sources in these fields all have either spectroscopic red-
shifts (spec-zs) or high-quality photo-zs, as have been exten-
sively examined in previous literature. Representative exam- 1 Defining ∆z = zphot − zspec , σNMAD is then the normalized median absolute
ples are listed in Column 7 of Table 1. Many more successful deviation of ∆z/(1 + zspec ), and outlier fraction is the fraction of sources
with |∆z| /(1 + zspec ) > 0.15. These two parameters are standard metrics
works built upon these redshifts have also indirectly justi- used to represent the photo-z quality.
fied their general reliability. When compared to the available
spec-zs, the photo-zs are of high quality – our sample has a
σNMAD of 0.03 (0.04) and an outlier fraction of 4% (15%) for
 5

in the hard band (HB)2 for CANDELS and the LSST DDFs. luminosity and LX ), BHAR can be expressed as follows.
The reason is to minimize the effects of obscuration. Select-
ing AGNs in soft bands (< 2 keV) is biased toward little or BHAR(M⋆ , z) =
Z +∞
no absorption. Since the obscuration level is known to be (1 − ϵ)kbol (M⋆ λ)M⋆ λ
p(λ | M⋆ , z)d log λ, (1)
correlated with λ (e.g., Ricci et al. 2017), soft-band-selected log λmin ϵc2
AGNs are expected to be biased in terms of λ. Besides, our
analyses need intrinsic LX , and HB fluxes are the least af- where ϵ is the radiative efficiency of the accretion. The key
fected by obscuration. To calculate LX and, consequently, step in measuring BHAR is hence to derive p(λ | M⋆ , z).
λ of these HB-detected sources, we use Equation A4 in Zou Some literature models the LX distribution instead of λ (e.g.,
et al. (2022) and adopt a photon index of 1.6. As discussed in Aird et al. 2012). These two approaches are equivalent, and
Yang et al. (2018), a photon index of 1.6 returns LX agreeing p(λ | M⋆ , z) and p(LX | M⋆ , z) are interchangeable. The
the best with those from X-ray spectral fitting. For eFEDS, as only reason for choosing one instead of the other is for con-
mentioned in Section 2.3, we use the de-absorbed 0.5−2 keV venience, as λ is a scaled parameter that can serve as a rough
flux in Liu et al. (2022) and convert it to LX assuming a pho- proxy for the Eddington ratio.
ton index of 1.8. Although eROSITA observations are more 3.1. Semiparametric Modeling of p(λ | M⋆ , z)
prone to obscuration effects, and it is less accurate to mea-
sure LX with soft X-rays below ≈ 2 keV, we have verified in We assume a double power-law with respect to λ for p(λ |
Appendix C that our median results remain similar when ex- M⋆ , z):
cluding eFEDS. It should be noted that we do not exclusively   −γ
λ
, λmin < λ ≤ λc
1
A λc


rely upon eFEDS to provide constraints at low-z and/or high- p(λ | M⋆ , z) = 

−γ (2)
A λ
 
, λ > λc
 2
M⋆ . The LSST DDFs, especially with the X-ray coverage 
λc
in Chen et al. (2018) and Ni et al. (2021b) added, already
have 12.6 deg2 of coverage with useful HB sensitivity (see The four parameters (A, λc , γ1 , γ2 ) are functions of (M⋆ , z).
Table 1), and thus can also provide beneficial constraints. We We require λc > λmin because, otherwise, the model will al-
define AGNs as those with log λ > log λmin = 31.5 and ne- ways degenerate to a single power law and has no depen-
glect the contribution of SMBHs with λ ≤ λmin to BHAR. dence on γ1 once λc lies below λmin . We also require γ2 > 0;
This is because few of the X-ray-detected AGNs are below otherwise, p(λ | M⋆ , z) will not be a probability measure,
λmin , and the emission from X-ray binaries may become non- and the model-predicted number of AGNs will diverge.3 It
negligible for low-λ sources. As we will show in Section 3.3, has been shown that a double power-law can indeed approxi-
BHAR is indeed dominated by sources above λmin . mate p(λ | M⋆ , z) well (e.g., Bongiorno et al. 2016; Aird et al.
In total, we have eight thousand AGNs and 1.3 million nor- 2018; Yang et al. 2018). Similarly, the observed AGN X-ray
mal galaxies, and they are plotted in the z − M⋆ and z − λ luminosity function (XLF) also follows a double power-law
planes in Figure 1, where each column presents fields with (e.g., Ueda et al. 2014), and a p(λ | M⋆ , z) roughly with
comparable depths and areas. Note that Yang et al. (2019), a double power-law shape is needed to reproduce the XLF
Zou et al. (2022), and Yu et al. (2023), from which our (Section 3.2).
adopted galaxy properties are taken, all have appropriately 3.1.1. The Detection Probability
considered the AGN emission for AGNs. We will also assess
the impact of AGNs that dominate the SEDs in Appendix D. We denote Pdet ( fX ) as the probability that a source with a
2 − 10 keV flux of fX is detected by a given X-ray survey.
Following Section 3.4 in Zou et al. (2023), we adopt the fol-
lowing functional form for Pdet ( fX ):
1
Pdet ( fX ) = erf b(log fX − a) + 1 ,


(3)
2
3. METHOD AND RESULTS
where a and b are parameters determining the shape of
Denoting p(λ | M⋆ , z) as the conditional probability den- Pdet ( fX ). We follow the same procedures as in Zou et al.
sity per unit log λ of a galaxy with (M⋆ , z) hosting an AGN (2023) to calibrate a and b and report the results in Ta-
with λ and kbol (LX ) as the LX -dependent 2 − 10 keV bolo- ble 1. Briefly, we compared the fX distribution with the
metric correction (i.e., the ratio between the AGN bolometric 2−10 keV log N −log S relation in Georgakakis et al. (2008),

2 The detection energy range for the HB has slightly different definitions in 3 Note that p(λ | M⋆ , z) is defined
R +∞ in the log λ space, and thus γ2 > 0 is
different fields – 2 − 7 keV for CANDELS and COSMOS, 2 − 12 keV for
W-CDF-S and ELAIS-S1, and 2 − 10 keV for XMM-LSS. sufficient and necessary for log λ p(λ | M⋆ , z)d log λ < +∞.
min
6

12.0 CANDELS LSST DDFs eFEDS

AGN
11.5
logM (M )

11.0
10.5
10.0
9.5
35.0
34.5
log (erg cm 2 M 1)

34.0
33.5
33.0
32.5
32.0
31.5
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
z z z
Figure 1. Our sample in the z − M⋆ (top) and z − λ (bottom) planes. The left, middle, and right panels are for CANDELS, the LSST DDFs, and
eFEDS, respectively. The points are AGNs. The background grey-scale cells in the left panel are the 2-D histogram of the number of normal
galaxies, with darker cells representing more galaxies. The apparent deficiency of sources in the high-z and/or low-M⋆ regime in the middle
and right panels is due to our M⋆ completeness cuts.

which is the well-determined expected surface number den- XLF-predicted intrinsic log N − log S relation is
sity per unit fX with the detection procedures deconvolved. Z +∞ Z 5
The comparison can return best-fit (a, b) parameters such that dVC
N( fX,int > S ) = A−1 all sky ϕL (LX , z) d log fX,int dz,
the convolution between the log N − log S relation and Pdet log S 0 dz
best matches the observed fX distribution. It is necessary to (4)
adopt a functional form because it improves the computa- fX,int
LX ( fX,int , z) = , (5)
tional speed by several orders of magnitude, as will be dis- η(z)
cussed below. The form of Equation 3 has been shown to be (1 + z)2−Γ
appropriate for X-ray surveys (e.g., Yan et al. 2023; Zou et al. η(z) = (6)
4πD2L
2023) because its overall shape is similar to X-ray sensitivity
curves, and, in our case, it indeed returns consistent BHAR where Aall sky is the all-sky solid angle, VC is the comoving
as in Yang et al. (2018), who did not adopt this functional volume within a redshift of z, η(z) is a function of z con-
form for Pdet . verting LX to the intrinsic 2 − 10 keV flux for a power-law
There is a subtle difference between eFEDS and the other X-ray spectrum with a power-law photon index of Γ = 1.8,
fields. For the latter, their fX is the observed value taken and DL is the luminosity distance. We limit the integration to
from the original X-ray catalogs. The log N − log S relation z < 5 because the contribution of higher-redshift sources to
is also for the observed fX ; thus, Pdet is for the observed fX . the total source number is negligible. Similarly, the predicted
For eFEDS, we adopt the intrinsic, de-absorbed 0.5 − 2 keV observed log N − log S relation is
flux in Liu et al. (2022) and multiply it by 1.57 to convert it Z +∞ Z 5 Z 24
to the intrinsic 2 − 10 keV flux assuming a photon-index of dVC
N( fX,obs > S ) = A−1 ϕL (LX , z)
Γ = 1.8. For consistency, we should correct the log N − log S all sky
log S 0 20 dz
relation such that it works for the intrinsic fX . We use the × p (NH | LX , z) d log fX,obs dzd log NH , (7)
XLF (ϕL ) in Ueda et al. (2014) to derive the correction. The
fX,obs
LX ( fX,obs , z, NH ) = , (8)
η(z)α(NH , z)
 7

where the NH function p (NH | LX , z) is the conditional prob- When compared with the observed data, the log-likelihood
ability density per unit log NH of an AGN with (LX , z), as function (e.g., Loredo 2004) is
given in Section 3 of Ueda et al. (2014). This function is nor-
R 24 ngal AGN
nX
malized such that 20 p (NH | LX , z) d log NH = 1. α(NH , z)
X
ln L = − T gal,s + ln p(λ s | M⋆,s , z s ), (9)
is the absorption factor for a source with Γ = 1.8 and is s=1 s=1
calculated based on photoelectric absorption and Compton- Z +∞
scattering losses (i.e., zphabs×cabs) in XSPEC. T gal = p(λ | M⋆ , z)Pdet ( fX (λM⋆ , z))d log λ, (10)
log λmin
The XLF-predicted N( fX,obs > S ) is similar to the
observed log N − log S relation, with an absolute dif- fX (λM⋆ , z) = λM⋆ η(z), (11)
ference generally below 0.2 dex. We found that
where η(z) is given in Equation 6. We adopt Γ = 1.8 and 1.6
log N( fX,int > S )/N( fX,obs > S ) is almost a constant around
 
for eFEDS and the other fields, respectively. Different Γ val-
0.15 dex at log S > 10−14 erg cm−2 s−1 , and thus we add
ues are adopted because the adopted fX inside our Pdet func-
0.15 dex to the observed log N−log S relation in Georgakakis
tion is the intrinsic value for eFEDS while being the observed
et al. (2008) to approximate the intrinsic relation. Applying
one for the other fields (Section 3.1.1). Equation 10 involves
this intrinsic relation for our calibration in Equation 3, we
an integration, and Equation 9 computes Equation 10 many
can obtain the eFEDS Pdet as a function of the intrinsic fX .
times in the summation for a single evaluation of L. Numer-
Given that the intrinsic fX instead of the observed fX is al-
ically integrating Equation 10 is slow, making it impractical
ways adopted in our analyses of eFEDS, the fact that eFEDS
to sample more than one or two dozen free parameters (as
is more easily affected by absorption has been appropriately
will be shown later, we will have 104 free parameters). For-
accounted for and absorbed into Pdet . For example, the fact
tunately, as previously suggested in Zou et al. (2023), Equa-
that obscured AGNs may be missed by eFEDS causes the a
tion 10 can be analytically solved when choosing appropriate
value to slightly shift to a larger value due to the correction
functional forms for p(λ | M⋆ , z) and Pdet ( fX ), and our Equa-
applied to the observed log N −log S relation. One may won-
tions 2 and 3 enable this. This is one of the most important
der why we convert the 0.5 − 2 keV flux to 2 − 10 keV flux
steps enabling our whole semiparametric analyses.
instead of directly using 0.5 − 2 keV flux. Since the intrinsic
We define
flux is always adopted, the conversion, in principle, would
not cause systematic biases. The main reason is that the cor- I(γ, λ1 , λ2 , A, λc ; M⋆ , z) =
rection to the log N − log S relation is considerably smaller Z log λ2 !−γ
λ
for the 2 − 10 keV band than for the 0.5 − 2 keV band. A Pdet (λM⋆ η(z))d log λ, (12)
log λ1 λc
One caveat is that we limit the integration range of NH in
Equation 7 to be below 1024 cm−2 , which equivalently means Using Equation 21 in Zou et al. (2023), Equation 12 can be
that we neglect the contribution from Compton-thick (CT) reduced as follows.
AGNs with NH > 1024 cm−2 for eFEDS. Similarly, in our !−γ !−γ
λ1 λ2
(
other fields observed by Chandra or XMM-Newton, we also A
I= [erf(x1 ) + 1] − [erf(x2 ) + 1]
implicitly neglect most CT AGNs because they can hardly 2γ ln 10 λc λc
−γ
be detected even in the HB. Generally, X-ray observations  10a− γ4bln210  "

 erf x1 + γ ln 10 − erf x2 + γ ln 10
! !# )
below 10 keV cannot provide effective constraints for the −  ,
λc M⋆ η  2b 2b
CT population, and the intrinsic fraction of CT AGNs is
highly uncertain (e.g., Ananna et al. 2019). Therefore, any (13)
xk = b log(λk M⋆ η) − a , k = 1, 2.
 
attempt to measure the intrinsic CT population properties us- (14)
ing X-rays below 10 keV is likely prone to large systematic
uncertainties. The CT population might indeed contribute to Equation 10 can then be expressed as follows.
the SMBH growth and is missed by our measurements, es-
T gal (A, λc , γ1 , γ2 ; M⋆ , z) = I(γ1 , λmin , λc , A, λc ; M⋆ , z)
pecially at high redshift (e.g., Yang et al. 2021), but obser-
vations in the regime insensitive to the CT obscuration are + I(γ2 , λc , +∞, A, λc ; M⋆ , z). (15)
necessary to reveal it (e.g., Yang et al. 2023).
The above equations express L as a function of
(A, λc , γ1 , γ2 ), which themselves are functions of (M⋆ , z).
The dependences of (A, λc , γ1 , γ2 ) on (M⋆ , z) lack clear
guidelines, and we use a nonparametric approach to model
them. We divide the (M⋆ , z) plane into N M × Nz grids and
3.1.2. The Likelihood adopt the (A, λc , γ1 , γ2 ) values in each grid element as free
parameters; i.e., we have 4N M Nz free parameters in total.
8

Such an approach is conceptually similar to and was indeed where X denotes each one of (log A, log λc , γ1 , γ2 ), and σX is
initially inspired by the “gold standard” nonparametric star- our chosen a priori parameters to quantify the overall varia-
formation history (e.g., Leja et al. 2019) in SED fitting. In tions of X across the whole fitting ranges. The goal of this
a strict statistical sense, a method is called “nonparametric” continuity prior is to transport information among grid ele-
only if the number of free parameters scales with the number ments. Without this prior, the fitted parameters in each grid
of data points. In contrast, we used a fixed number of free element become unstable and vary strongly. This prior is de-
parameters, which does not exactly satisfy the statistical def- fined in a way such that the information flow is roughly inde-
inition. Although we can easily adjust N M and Nz so that the pendent of the grid size. The continuity prior is defined only
number of free parameters scales with the number of data for the differences, and we need a further prior for the X’s in
points, this makes the computation infeasible because we a single cell and adopt it as flat in the (log A, log λc , γ1 , γ2 )
have millions of galaxies; besides, with our continuity prior space. We set bounds for these parameters to ensure pro-
in Section 3.1.3, further increasing the number of free pa- priety of the prior (Tak et al. 2018): −10 < log A < 10,
rameters does not improve our results materially. In our con- log λmin < log λc < 40, −5 < γ1 < 10, and 0 < γ2 < 10.
text, we use the word “nonparametric” because our number These ranges are sufficiently large to encompass any reason-
of free parameters is far larger than that of typical parametric able parameter values. Our posterior (Section 3.1.4) may also
methods, and our method is effectively similar to the fully become less numerically stable outside these bounds. The re-
nonparametric approach. This same argument also works sulting prior is explicitly shown below.
for “nonparametric star-formation history” in, e.g., Leja et al.  Nz
NX M −1 X
(2019). 1 X  (Xi+1, j − Xi j )2
ln πcont = − N M
This method has an important advantage over a parametric 2 X i=1 j=1
σ2X
one in our case. As Figure 1 shows, most of our data are clus- NM N z −1

tered within a small region of the (M⋆ , z) plane – the number
X X (Xi, j+1 − Xi j )2 
+Nz  . (19)
of sources significantly decreases at both low z (≲ 0.8) and i=1 j=1
σ2X
high z (≳ 2), the number of galaxies strongly depends on
Note that it is defined in the (log A, log λc , γ1 , γ2 ) space, and
M⋆ , and most AGNs are confined within 1010.5 ≲ M⋆ ≲
an appropriate Jacobian determinant should be added when
1011.2 M⊙ . This indicates that if we assume any parametric
transforming the parameter space. For sampling purposes,
form for (A, λc , γ1 , γ2 ), the fitted parameters will be domi-
variable transformations are usually needed.
nated by the small but well-populated region in the (M⋆ , z)
We rely on previous literature to set appropriate values for
plane. Especially, one strong argument disfavoring paramet-
σX . Yang et al. (2018) used a double power-law similar to
ric fitting is that our ultimate goal is to derive BHAR across
ours to fit p(λ | M⋆ , z), and their best-fit parameters (see their
all redshifts, but any parametric fitting will return results
Equation 16) span ranges of −3.53 < log A < −0.86, 31.73 <
dominated by sources in a small redshift range (e.g., Yang
log λc < 34.98, γ1 = 0.43, 1.55 < γ2 < 3.55 across our pa-
et al. 2018). Our semiparametric settings avoid this problem.
rameter spaces. Bongiorno et al. (2016) modeled the bivari-
Equation 9 then becomes
ate distribution function of M⋆ and λ for AGNs, which can be
Nz
NM X
X h gal converted to p(λ | M⋆ , z) by dividing it by the galaxy stellar
ln L = −ni j T gal (Ai j , λc,i j , γ1,i j , γ2,i j ; M⋆,i , z j )
mass function (SMF), and the corresponding p(λ | M⋆ , z) is
i=1 j=1
also a double power-law. We use the SMF in Wright et al.
nAGN
Xij
i (2018) for the conversion, and the best-fit double power-
+ ln pi j (λ s | M⋆,s , z s ) , (16) law parameters in Bongiorno et al. (2016) span ranges of
s=1
−5.28 < log A < −0.08, 33.32 < log λc < 34.52, −0.67 <
gal
where ni j and nAGNij are the numbers of galaxies and AGNs γ1 < 1.62, γ2 = 3.72. Aird et al. (2018) nonparametrically
within the (i, j) bin, respectively. L is defined for each indi- modeled p(λ | M⋆ , z), and we use our double power-law
vidual survey field, and they are added together to return the model to fit their results above M⋆ = 109.5 M⊙ by mini-
final likelihood. mizing the Kullback-Leibler divergence of our model rela-
tive to theirs. The returned best-fit values range between
3.1.3. The Prior
−2.87 < log A < −0.69, 31.84 < log λc < 34.04, −0.58 <
We adopt a continuity prior: γ1 < 0.52, 0.72 < γ2 < 1.67. Another independent way to
 σ2X 
  estimate p(λ | M⋆ , z) is based on the fact that p(λ | M⋆ , z),
Xi+1, j − Xi j ∼ N 0,  , (17) by definition, can predict the XLF when combined with the
NM
SMF (see Equation 22 and Section 4.1 for more details). We
 σ2 
 
Xi, j+1 − Xi j ∼ N 0, X  , (18) estimate parameters of p(λ | M⋆ , z) such that, when using the
Nz SMF in Wright et al. (2018), the predicted XLF can match
 9

the best with the XLF in Ueda et al. (2014). This returns that the true dependence is indeed well-approximated by a
−2.81 < log A < −1.08, 32.72 < log λc < 33.77, −0.35 < double power-law when λ > λmin . Previous works have
γ1 < 0.90, 2.46 < γ2 < 2.82. Taking the union of these shown that a double power-law indeed works, and, as far as
estimations, the ranges should span no more than −5.28 < we know, there is no clear evidence suggesting that this as-
log A < −0.08, 31.73 < log λc < 34.98, −0.67 < γ1 < sumption would fail. Especially, the nonparametric form of
1.62, 0.72 < γ2 < 3.72. We adopt σX as one third of the p(λ | M⋆ , z) inferred from Aird et al. (2018) is also roughly
widths,4 i.e., σlog A = 1.7, σlog λc = 1.1, σγ1 = 0.8, σγ2 = 1.0. a double power-law. The adopted approach essentially de-
In fact, our prior setting is essentially a rasterized approxi- pends on our ultimate goal. It is certainly better to mini-
mation to the continuous surface of a Gaussian process (GP) mize the assumption for p(λ | M⋆ , z) and adopt a nonpara-
regression (e.g., Rasmussen & Williams 2006). This is be- metric form for it if the ultimate goal is to derive the shape
cause the blocky prior surface over the (M⋆ , z) plane becomes of p(λ | M⋆ , z). However, our goal is different – we are ulti-
the non-parametric GP-based continuous surface as the res- mately interested in BHAR and thus want to assume a double
olution of the grid increases (i.e., increasing N M and Nz to power-law form for p(λ | M⋆ , z).
the infinity). Therefore, a full GP regression involves a large
number of free parameters scaling with the galaxy sample 3.2. Hamiltonian Monte Carlo Sampling of p(λ | M⋆ , z)
size (≈ 106 ), while our rasterized approach only involves 104 Given the high dimensionality, Hamiltonian Monte Carlo
parameters. GP also requires computations of O(n3 ) for ma- (HMC; e.g., Betancourt 2017) should be one of the most
trix inversions, while our approach turns the matrix-inversion practical methods to sample the posterior. As far as we
problem into products of multiple univariate Gaussian densi- know, other sampling methods either cannot work efficiently
ties. Due to these reasons, a full GP regression is computa- in our high-dimension case (e.g., the traditional Metropo-
tionally infeasible in our case, but our approach effectively lis–Hastings algorithm) or do not have well-developed pack-
works similarly and is much less computationally demand- ages readily available (e.g., Bayer et al. 2023). HMC needs
ing. both the posterior and its gradient in the parameter space.
The posterior has been presented in the previous subsec-
3.1.4. The Posterior tions, and we present the gradient in Appendix A. We use
The posterior is DynamicHMC.jl5 to conduct the HMC sampling. We adopt
X N M = 49 and Nz = 50. The sampling results are presented in
ln P = ln L + ln πcont . (20) Figure 2. These parameter maps will be released online.
field To examine our fitting quality, we compare the model
We call our overall modeling “semiparametric” because we p(λ | M⋆ , z) with the observed values. We use the nobs /nmdl
adopt p(λ | M⋆ , z) as a parametric function of λ, while the method to obtain binned estimators of p(λ | M⋆ , z), as out-
dependences of (A, λc , γ1 , γ2 ) on (M⋆ , z) are nonparametric. lined in Aird et al. (2012). For a given (z, M⋆ , λ) bin ranging
Readers may wonder why we do not also adopt a nonpara- [zlow , zhigh ] × [M⋆,low , M⋆,high ] × [λlow , λhigh ], we denote the
metric function for p(λ | M⋆ , z). In principle, it could be done number of observed AGNs as nobs and calculate the model-
and was presented in Georgakakis et al. (2017) and Aird et al. predicted number as nmdl :
(2018). Since any model contains subjective assumptions, X Z log λhigh
the choice of the methodology should be guided by the as- nmdl = p(λ | M⋆,s , z s )Pdet ( fX (λM⋆,s , z s ))d log λ,
log λlow
sumptions we want to retain or avoid. Compared to nonpara- s

metric modeling, the assumptions of parametric models are (21)

much stronger. We nonparametrically model (A, λc , γ1 , γ2 )
where the summation runs over all the sources within
as functions of (M⋆ , z) because we genuinely do not know
[zlow , zhigh ]×[M⋆,low , M⋆,high ]. The observed binned estimator
their dependencies and thus want to minimize assumptions.
of p(λ | M⋆ , z) is then the fitted model evaluated at the bin
However, we are satisfied with and thus want to retain the
center scaled by nobs /nmdl , and its uncertainty is calculated
inherent assumption of our parametrization of p(λ | M⋆ , z)
from the Poisson error of nobs following Gehrels (1986). We
present our model p(λ | M⋆ , z) and the binned estimators
4 A nominal σ is often approximated by 1/4 of the range, according to the so- in Figure 3, and they are consistent. The uncertainties be-
called range rule of thumb. We have two dimensions
 √  in our case, and thus come large especially in the high-z/low-M⋆ and low-z/high-
the one-dimension σ can be chosen as 1/ 4 2 of the range. However, M⋆ regimes because of a limited number of AGNs being
we would like to be slightly more conservative. The reason is that previous
works mostly do not cover a parameter space as large as this work, and thus available. In the high-z/low-M⋆ regime, most of the con-
extrapolations are employed when computing the ranges. Some conserva-
tiveness can enable more flexibility to accommodate possible systematic
extrapolation errors in regimes not well-covered by previous works. 5 https://www.tamaspapp.eu/DynamicHMC.jl/stable/
 10
Median logA Median log c Median 1 Median 2
12.0 12.0 12.0 12.0
1.0 33.6 0.6 2.1
11.5 11.5 11.5 11.5
33.4 2.0
logM (M )

logM (M )

logM (M )

logM (M )
11.0 1.5 11.0 11.0 0.5 11.0
33.2 1.9
10.5 10.5 33.0 10.5 0.4 10.5 1.8
2.0
10.0 10.0 32.8 10.0 10.0 1.7
2.5 0.3
9.5 9.5 32.6 9.5 9.5 1.6
0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.2 0.3 0.4 0.5 0.6
log(1 + z) log(1 + z) log(1 + z) log(1 + z)
Uncertainty of logA Uncertainty of log c Uncertainty of 1 Uncertainty of 2
12.0 0.225 12.0 12.0 12.0
0.16 0.12
11.5 0.200 11.5 11.5 11.5 0.16
0.14 0.11
0.175
logM (M )

logM (M )

logM (M )

logM (M )
11.0 11.0 11.0 0.10 11.0
0.12 0.14
0.150
10.5 10.5 0.10 10.5 0.09 10.5
0.125 0.12
10.0 10.0 0.08 10.0 0.08 10.0
0.100
9.5 9.5 9.5 0.07 9.5 0.10
0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.2 0.3 0.4 0.5 0.6
log(1 + z) log(1 + z) log(1 + z) log(1 + z)

Figure 2. The sampled maps of (A, λc , γ1 , γ2 ). The top panels are the median posteriors, and the bottom panels are the 1σ uncertainties, defined
as the half width of the posterior’s 16th-84th percentile range.

straints are from deep CANDELS fields, especially GOODS- As a simple check, for the parameters in Figure 2, if we ex-
S, because the other fields are not sufficiently deep in both trapolate the integration in Equation 22 to (−∞, +∞), the typ-
X-rays and other multi-wavelength bands. For example, 60% ical ϕL,mdl will only increase by 0.01 dex at 43 < LX < 43.5,
(80%) of AGNs in our sample with M⋆ < 1010 M⊙ and z > 2 the lowest LX bin that we will later adopt in Section 4.1. This
(z > 3) are from GOODS-S. At z < 0.5, ≳ 60% of AGNs are increment is even smaller for higher LX bins.
from eFEDS, and even the 60 deg2 eFEDS is not sufficiently
large to effectively sample high-M⋆ sources at low redshift. 3.3. Measuring BHAR
We also plot several p(λ | M⋆ , z) results from previous works
Equation 1 converts p(λ | M⋆ , z) to BHAR. We adopt
and leave more detailed discussions on the comparison be-
ϵ = 0.1 and kbol from Equation 2 in Duras et al. (2020). In
tween our p(λ | M⋆ , z) and previous ones to Section 4.2.
principle, ϵ may depend upon other factors such as the accre-
As another independent check, p(λ | M⋆ , z), by definition,
tion state (e.g., Yuan & Narayan 2014), but it is infeasible to
can connect the SMF (ϕ M ) and XLF (ϕL ). That is, the SMF
accurately measure ϵ for our individual sources. We adopt ϵ
and p(λ | M⋆ , z) can jointly predict the XLF (e.g., Bongiorno
as 0.1 because it is a typical value for the general AGN pop-
et al. 2016; Georgakakis et al. 2017):
ulation (e.g., Brandt & Alexander 2015) and has been widely
Z log M⋆,max
used in previous literature (e.g., Yang et al. 2017, 2018, 2019;
ϕL,mdl (LX , z) = p(λ | M⋆ , z)ϕ M d log M⋆
log M⋆,min Ni et al. 2019, 2021a). The kbol relation in Duras et al. (2020)
Z log M⋆,max diverges at high LX . To avoid it, we cap kbol not to exceed
= p(LX /M⋆ | M⋆ , z)ϕ M d log M⋆ . 363, the value when the bolometric luminosity is 1014.5 L⊙ ,
log M⋆,min
which is roughly the brightest sample used in Duras et al.
(22)
(2020) to calibrate the kbol relation. We show the LX − kbol
Comparing ϕL,mdl and the observed XLF, ϕL,obs , can thus fur- relation in Figure 5, in which we also plot the relation used
ther assess our fitting quality. We adopt ϕ M in Wright et al. in Yang et al. (2018), derived from Lusso et al. (2012), for a
(2018) and the median parameter maps in Figure 2 to calcu- comparison. The two relations are similar, with a small off-
late ϕL,mdl . We present the comparison between ϕL,mdl and set of ≈ 0.07 dex that is almost negligible compared to the
ϕL,obs from Ueda et al. (2014) in Figure 4, and they agree BHAR uncertainty (Figure 6). The deviation of the two rela-
well. Note that for comparison purposes here, we do not tions at LX ≳ 1045 erg s−1 has little impact on BHAR because
need to optimize the computation of Equation 22; however, BHAR has little contribution from log λ ≳ 35 (see Figure 3).
we will present a more optimized computation algorithm Equation 1 ignores the contribution to BHAR from sources
later in Section 4.1, where we do need fast computational at log λ < 31.5 because X-ray binaries may not be negligible
speed. Also, note that Equation 22 ignores the contribution at lower λ, and our X-ray surveys can hardly provide strong
from sources with M⋆ below M⋆,min = 109.5 M⊙ or above constraints in the low-λ regime. However, this will not cause
M⋆,max = 1012 M⊙ to the XLF. This is appropriate because material biases because BHAR is dominated by sources at
the XLF is dominated by AGNs with 109.5 < M⋆ < 1012 M⊙ . log λ ≳ 31.5 (e.g., Section 3.2.3 in Yang et al. 2018). We
 11

Aird+12 Aird+18 Bongiorno+16 Georgakakis+17 Yang+18 This work

9.5 < logM < 10.0 10.0 < logM < 10.5 10.5 < logM < 11.0 11.0 < logM < 11.5 11.5 < logM < 12.0
0
logp( M , z)

0.1 < z < 0.5

2
4
0
logp( M , z)

0.5 < z < 1.0

2
4
0
logp( M , z)

1.0 < z < 1.5

2
4
0
logp( M , z)

1.5 < z < 2.0

2
4
0
logp( M , z)

2.0 < z < 2.5

2
4
0
logp( M , z)

2.5 < z < 3.0

2
4
0
logp( M , z)

3.0 < z < 4.0

2
4
32 33 34 35 32 33 34 35 32 33 34 35 32 33 34 35 32 33 34 35
log log log log log
Figure 3. Comparison between our p(λ | M⋆ , z) and other measurements. The red points are the binned estimators with 1σ error bars based
on our data. The blue curves are our fitted median p(λ | M⋆ , z), evaluated at the bin centers, and the blue shaded regions are the corresponding
90% confidence ranges. The black solid straight lines are the single power-law models in Aird et al. (2012). The dash-dotted and dashed curves
are the double power-law models in Bongiorno et al. (2016) and Yang et al. (2018), respectively. The cyan and orange shaded regions are the
90% confidence intervals of the nonparametric p(λ | M⋆ , z) in Aird et al. (2018) and Georgakakis et al. (2017), respectively.
 12

2 z = 0.1 z = 0.3 z = 0.5 z = 0.7

log L (Mpc 3 dex 1) log L (Mpc 3 dex 1) log L (Mpc 3 dex 1)

4
6
8
L, mdl
10 Ueda et al. (2014)

2 z = 0.9 z = 1.1 z = 1.4 z = 1.8

4
6
8
10
2 z = 2.2 z = 2.7 z = 3.5
4
6
8
10
42 43 44 45 46 42 43 44 45 46 42 43 44 45 46
logLX (erg s 1) logLX (erg s 1) logLX (erg s 1)
Figure 4. The XLFs at different redshifts. The red (blue) data points are the soft-band (hard-band) XLFs in Ueda et al. (2014). The cyan curves
are the best-fit XLF models in Ueda et al. (2014), and the black curves are our ϕL,mdl based on the median parameter maps in Figure 2 and the
SMF in Wright et al. (2018). The absorption correction has been appropriately applied for both our measurements (see Section 3.1.1) and the
XLFs in Ueda et al. (2014). Our models agree with the observed XLFs well.
 13

This work 12.0 Median logBHAR

2.5 Yang et al. (2018) 0
11.5
2.0 1

logM (M )
11.0
logkbol

2
1.5
10.5 3
1.0 10.0 4

0.5 9.5 5
42.0 42.5 43.0 43.5 44.0 44.5 45.0 45.5 46.0 0.1 0.2 0.3 0.4 0.5 0.6
logLX (erg s 1) log(1 + z)

Figure 5. The adopted LX − kbol relation, taken from Duras et al. 12.0 Uncertainty of logBHAR
(2020). The adopted relation used in Yang et al. (2018), which is 0.30
adjusted from Lusso et al. (2012), is also plotted for comparison. 11.5
0.25

logM (M )
have also tried pushing the lower integration bound in Equa- 11.0 0.20
tion 1 down by 2 dex, and the returned BHAR would only
increase by a typical value of ≈ 0.02 dex and no more than 10.5 0.15
≈ 0.1 dex. Such a difference is much smaller than the fit- 0.10
ted BHAR uncertainty. This exercise may even overesti- 10.0
mate the influence because p(λ | M⋆ , z) may bend downward 0.05
or quickly vanish at very small λ (Aird et al. 2017, 2018). 9.5 0.00
Therefore, the cut at log λ = 31.5 is not expected to cause 0.1 0.2 0.3 0.4 0.5 0.6
material biases to BHAR. log(1 + z)
We show our sampled BHAR results in Figure 6, and the
BHAR maps will be released online. The median map clearly 1 z = 0.1 z = 1.5
z = 0.5 z = 2.0
shows that BHAR increases with both M⋆ and z, qualitatively
0 z = 0.8 z = 3.0
logBHAR (M yr 1)

consistent with the conclusions in Yang et al. (2018). The un- z = 1.0 z = 4.0
certainty map reveals that the BHAR constraints at both the 1
low-z/high-M⋆ and the high-z/low-M⋆ regime are relatively
more limited. We will present more quantitative comparison 2
with Yang et al. (2018) and other works in Section 4.2. Be-
sides, we verified that AGN-dominated sources do not cause 3
material biases to our BHAR measurements in Appendix D. 4
There are slight, local fluctuations in BHAR that are
caused by the statistical noise of the data and are confined 5
within the extent allowed by our prior, and the BHAR map is 9.5 10.0 10.5 11.0 11.5 12.0
smooth globally, as can be seen in the top panel of Figure 6.
logM (M )
The fluctuation levels and BHAR uncertainties depend on our
Figure 6. The top and middle panels are the sampled BHAR maps,
prior settings but almost not on our bin size because our bins where the unit of BHAR is M⊙ yr−1 . The bottom panel shows the
are set to be correlated (Section 3.1.3). For example, relax- BHAR − M⋆ relation at several redshifts, where the solid curves are
ing the prior by choosing larger σX would return larger fluc- the median values, and the shaded regions are the corresponding 1σ
tuations and uncertainties. This “arbitrariness” is inherent in uncertainty ranges. BHAR generally increases with both M⋆ and z.
modeling.6 Overall, our prior is reasonably constructed (Sec-
tion 3.1.3) and provides beneficial regularizations. We have
6 For the widely used method of binning the parameter space and assuming
assessed the potential issue of whether such arbitrary choices
each bin is independent, there is a similar arbitrariness in choosing the bin affect the following posterior inferences and the resulting sci-
size, and the uncertainties in this case would depend on the bin size. entific conclusions qualitatively. For example, we have con-
14

ducted a sensitivity check of our priors and confirmed that to the integration in Equation 22 is
the impact of lower or higher resolution of the prior surface
(corresponding to larger or smaller bin sizes) does not influ- ψDP (A, λc , γ1 , γ2 , M1 , M2 ; LX )
Z log M2
ence the resulting posterior inference in a noticeable way, and
= p(LX /M⋆ | M⋆ , z)ϕ M d log M⋆
changing σX generally would not cause material changes of log M1
the median BHAR map. 




ψ(γ2 , M1 , M2 , A, λc ; LX ), λc ≤ LX /M2

ψ(γ2 , M1 , LX /λc , A, λc ; LX )



=

4. DISCUSSION
+ψ(γ1 , LX /λc , M2 , A, λc ; LX ), LX /M2 < λc < LX /M1 ,




Given that this article is already lengthy and full of techni-



ψ(γ , M , M , A, λ ; L ), λ ≥ L /M


cal details, we decide to present more scientific investigations 1 1 2 c X c X 1

of our results in future dedicated works. However, we would (25)

like to present brief, immediate, but sufficiently informative
Equation 22 is thus
explorations of our results in this section, which helps jus-
tify the quality and serves as a precursor of further detailed NM
X
scientific investigations. ϕL,mdl = ψDP (Ai jz , λc,i jz , γ1,i jz , γ2,i jz , MLB,i , MLB,i+1 ; LX ),
i=1
(26)
4.1. Adding External Constraints from the SMF and XLF where log MLB,i = log M⋆,min + (i − 1)/N M ×
Section 3.2 uses the SMF and XLF to examine the fitting log(M⋆,max /M⋆,min ) is the lower bound of the ith M⋆ -grid
quality of p(λ | M⋆ , z). It is also possible to follow a reversed element, and jz is the index of the z-grid element containing
direction – we can add external constraints from the SMF and z.
XLF (named “the SMF-XLF constraints” hereafter) into our We then follow the procedure in Yang et al. (2018) to com-
posterior. This approach was adopted in Yang et al. (2018). pare ϕL,mdl and ϕL,obs in Ueda et al. (2014). ϕL,obs is eval-
As a start, we revisit the numerical computations of ϕL,mdl uated at several (LX , z) values, and the number of sources
in Equation 22. Again, numerical integrations should be (nXLF ) in Ueda et al. (2014) contributing to ϕL,obs is recorded.
avoided whenever possible, and we hence derive an analyti- Following Yang et al. (2018), we use the soft-band XLF at
cal formula for ϕL,mdl . ϕ M is expressed as a two-component LX > 1043 erg s−1 in Ueda et al. (2014). Their soft-band
Schechter function in Wright et al. (2018): XLF has been corrected for obscuration and spans a wider
LX range compared to their HB XLF, and their soft-band and
 !α1 +1 !α +1  HB XLFs are also consistent (see Figure 4). The LX cut at
dn M⋆
ϕ1 M⋆ M 2 
ϕM = = ln 10 e− Mc + ϕ2  ,

1043 erg s−1 is adopted to avoid contamination from X-ray bi-
d log M⋆ Mc Mc
naries. The log-likelihood when comparing ϕL,mdl and ϕL,obs
(23) is
ϕL,mdl,k XLF
X ! !
where (Mc , α1 , α2 , ϕ1 , ϕ2 ) are redshift-dependent parameters. ln LSMF−XLF = ln Pr Poisson nk = nXLF
ϕL,obs,k k
We further define an auxiliary function ψ such that the model- k
ϕL,mdl,k ϕL,mdl,k
" ! #
predicted XLF in Equation 22 can be simplified as summa-
X
= nXLF ln − + const.,
tions of ψ (see below). k
k
ϕL,obs,k ϕL,obs,k
(27)
ψ(γ, M1 , M2 , A, λc ; LX ) ϕL,mdl,k = ϕL,mdl (LX,k , zk ), (28)
Z log M2 !−γ
LX where k runs over all the LX and z bins of the observed XLF
= A ϕ M d log M⋆
log M1 M⋆ λc in Ueda et al. (2014). This term is called the SMF-XLF like-
!−γ "
lihood in Yang et al. (2018).
!
LX M1 M2
=A ϕ1 ΓGI α1 + γ + 1, ,
M⋆ λc Mc Mc To add the SMF-XLF constraints, Equation 20 should be
M1 M2
!# modified as follows.
+ϕ2 ΓGI α2 + γ + 1, , , (24) X
Mc Mc ln P = ln L + ln LSMF−XLF + ln πcont . (29)
field
x2 ζ−1 −t
where ΓGI (ζ, x1 , x2 ) =
R
x1
is the generalized incom-
t e dt Its gradient is presented in Appendix B for HMC sampling.
plete Gamma function. The contribution of each grid element We then sample the above posterior with HMC and present
 15

the resulting BHAR in Figure 7. The BHAR curves with p(λ | M⋆ , z) differences generally have limited data and are
or without the SMF-XLF constraints are largely consistent far away from the bulk of other data points, and the results
with a small (< 1σ) difference. This is expected because in these regions are subject to large uncertainties. For Bon-
Figure 4 shows that our BHAR without the SMF-XLF con- giorno et al. (2016), their p(λ | M⋆ , z) is similar to ours at
straints leads to consistent XLFs with those in Ueda et al. 10 ≲ log M⋆ ≲ 11.5 and 1.5 ≲ z ≲ 2 but has a much steeper
(2014). low-λ power-law index at z < 1.5. Two reasons may be
Although there is a good consistency after adding the responsible for the difference – the data used in Bongiorno
SMF-XLF constraints in our case, extra cautions are gen- et al. (2016) are not sufficiently deep to effectively probe
erally needed. The SMF and XLF taken from other litera- the low-λ regime; their model always fixes the breakpoint
ture works usually involve inherent assumptions about their at log λ = 33.8 when M⋆ = 1011 M⊙ , while our results sug-
parametric forms. When putting the SMF and XLF into our gest that the breakpoint tends to become smaller as redshift
posterior, we will inevitably “absorb” these assumptions. Be- decreases.
sides, the original data used to measure the XLF may over- Georgakakis et al. (2017) and Aird et al. (2018) adopted
lap with one’s dataset, especially given that the X-ray data in nonparametric methods to measure p(λ | M⋆ , z) without as-
GOODS-S are also necessary to constrain the XLF at low-LX suming a double power-law form. Our results show good
and/or high-z. Such an overlap causes double counting of the agreement with theirs, especially in regimes well-covered by
involved sources. Especially, more considerations would be the data, suggesting that a double power-law is indeed a good
needed if the posterior is dominated by the SMF-XLF con- approximation of p(λ | M⋆ , z). Nonetheless, some differ-
straints. ences are worth noting. At log λ ≳ 34 where the data be-
come limited, the p(λ | M⋆ , z) in Aird et al. (2018) tends to
4.2. Comparison with Previous Works be flatter than ours, while that in Georgakakis et al. (2017)
tends to be steeper than ours. This high-λ regime is highly
Figure 3 compares our p(λ | M⋆ , z) with some repre-
uncertain and subject to the adopted methodology – for in-
sentative results in previous literature. Aird et al. (2012;
stance, the prior adopted in Aird et al. (2018) prefers a flat-
black solid lines in Figure 3) used a single power-law to
ter slope at high λ. Another feature is that the p(λ | M⋆ , z)
fit p(LX | M⋆ , z) at z < 1, which is converted to a sin-
in Aird et al. (2018) sometimes shows downward bending
gle power-law p(λ | M⋆ , z) in Figure 3. The single power-
at log λ ≈ 32 − 33, while that in Georgakakis et al. (2017)
law curves broadly follow our double power-law ones, and
does not show a clear bending, although the large uncertainty
the single power-law index lies within the range between γ1
may be responsible for the lack of bending. In principle,
and γ2 . This indicates that a single power-law model can
a downward bending at some low λ is expected; otherwise,
serve well as the first-order approximation of p(λ | M⋆ , z),
p(λ | M⋆ , z) would diverge. Such bending can also be seen
as has been widely adopted in other works (e.g., Bongiorno
in Georgakakis et al. (2017), but below log λmin = 31.5 (see,
et al. 2012; Wang et al. 2017; Birchall et al. 2020, 2023;
e.g., their Figure 7). Our double power-law model is unable
Zou et al. 2023), especially when the data are limited. How-
to capture this feature, and Figure 3 shows that the bending in
ever, the real p(λ | M⋆ , z) is more complicated, and a double
Aird et al. (2018) mainly becomes prominent at high redshift
power-law model can return better characterizations. As Fig-
(z ≳ 3).
ure 3 shows, the binned p(λ | M⋆ , z) estimators generally do
Another metric that can be measured from p(λ | M⋆ , z)
not show systematic deviations from our double power-law
is the fraction of galaxies hosting accreting SMBHs above a
curves (e.g., no further breaks are visible), and thus a double
given λ threshold (λthres ), as calculated below.
power-law model is sufficient to capture the main structures
of p(λ | M⋆ , z) at λ > λmin .
Z +∞
Bongiorno et al. (2016) and Yang et al. (2018) adopted a fAGN (λ > λthres ) = p(λ | M⋆ , z)d log λ. (30)
log λthres
double power-law model similar to ours, and we plot their
results as the dash-dotted and dashed lines in Figure 3, re- For a consistent comparison with Aird et al. (2018), we adopt
spectively. Our p(λ | M⋆ , z) curves are nearly identical to the same λthres = 32 as theirs. We calculate fAGN at several
those in Yang et al. (2018) at 10 ≲ log M⋆ ≲ 11.5 and (M⋆ , z) values and plot the results in Figure 8. Our results
1 ≲ z ≲ 2.5 but begin diverging in other parameter ranges. In generally agree well with those in Aird et al. (2018) and fol-
the lowest-mass bin (9.5 < log M⋆ < 10), our p(λ | M⋆ , z) is low similar evolutionary trends with respect to M⋆ and z. At
still similar to those in Yang et al. (2018) at log λ ≲ 33.5 log M⋆ ≳ 10, fAGN increases with z up to z ≈ 1.5 − 2 and
but is lower than theirs at higher λ. In the highest-mass reaches a plateau at higher redshift; while for less-massive
bin (11.5 < log M⋆ < 12), our p(λ | M⋆ , z) is larger at galaxies, the redshift evolution is weaker. At low redshift
z ≲ 2 and smaller at z ≳ 2 than for Yang et al. (2018). It (z ≲ 0.5), fAGN is similar regardless of M⋆ , and this con-
should be noted that these parameter regions with noticeable clusion can be further extended down to log M⋆ < 9.5, as
16

1
z = 0.1 z = 0.3 z = 0.5
0
logBHAR (M yr 1)

1
2
3
4
5
1
z = 0.8 z = 1.0 z = 1.5
0
logBHAR (M yr 1)

1
2
3
4
5
1
z = 2.0 z = 3.0 z = 4.0
0
logBHAR (M yr 1)

1
2
3
4
5
9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0
logM (M ) logM (M ) logM (M )
Figure 7. BHAR as a function of M⋆ at several redshifts. The red curves are our median BHAR with the SMF-XLF constraints added, and the
orange shaded regions represent the corresponding 1σ and 2σ uncertainty ranges. The blue curves are our median BHAR and 1σ uncertainties
without the SMF-XLF constraints.
Zou et al. (2023) showed that the λ-based fAGN in the dwarf- to theirs, but some subtle differences exist. Our low-mass
galaxy population in this redshift range is also similar to BHAR at log M⋆ ≲ 10 is slightly smaller across all red-
fAGN in massive galaxies. At higher redshift (z ≳ 1), the shifts, though not very significant. Our high-mass BHAR at
dependence of fAGN on M⋆ becomes more apparent due to log M⋆ ≳ 11.5 differs the most with Yang et al. (2018), and
M⋆ -dependent redshift evolution rates of fAGN , and there is a ours tends to be smaller at z ≳ 3 while being larger at z ≲ 2.
positive correlation between fAGN and M⋆ at log M⋆ ≲ 10.5. These differences originate from different p(λ | M⋆ , z), as
However, for massive galaxies with log M⋆ ≳ 10.5, fAGN discussed earlier in this section. As shown in Figure 3, our
nearly does not depend on M⋆ . A full physical explanation low-mass p(λ | M⋆ , z) is smaller than for Yang et al. (2018)
of these complicated correlations between fAGN and (M⋆ , z) only at high λ, and thus the low-mass BHAR difference is
will require further detailed analyses of p(λ | M⋆ , z) with at small. Our high-mass p(λ | M⋆ , z) at log M⋆ ≳ 11.5, instead,
least partially physically driven modeling, and we leave these shows a redshift-dependent difference in the normalization.
analyses for future work. Nevertheless, the uncertainties in these extreme regimes are
We further compare our BHAR with those from Yang et al. large, and they are also subject to model choices. Dedicated
(2018) in Figure 9. Our median relation is largely similar analyses of these extreme-mass sources with deeper or wider
 17

star formation and SMBH accretion. The BHAR difference

between the two populations also depends on M⋆ and z. At
log M⋆ ≲ 10.5, the difference is generally small across most
10 1 of the redshift range. At higher mass, the difference is small
at low redshift but becomes apparent when z increases to 1
and further decreases at higher redshift. There is also tenta-
fAGN

tive evidence suggesting that the redshift at which the differ-

ence reaches its peak might also shift with M⋆ , with the peak
of the BHAR difference of higher-mass galaxies occurring at
10 2
logM = 9.75 logM = 11.25 higher redshift.
logM = 10.25 logM = 11.75 One caveat that should be noted is that our results de-
logM = 10.75 pend on the classification between star-forming and quies-
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 cent galaxies. Such a classification is more reliable at M⋆ ≲
z 1010.5 M⊙ , but it may become sensitive to the adopted method
at higher M⋆ and lower z (e.g., Donnari et al. 2019). Cristello
Figure 8. fAGN evaluated at several (M⋆ , z) values versus z. Our et al. (2024) show that the star-forming and quiescent sub-
results are plotted as open circles with 1σ error bars, where different populations cannot be safely defined for massive galaxies,
colors represent different M⋆ . The dashed lines are those in Aird
and Feldmann (2017) also argued that the bimodal separa-
et al. (2018).
tion is not necessarily appropriate. The proposed redshift-
dependent maximum M⋆ values for reliable classifications in
data may be necessary to further pin down the uncertainty.
Cristello et al. (2024) can be well described by the following
Another important difference is that the BHAR − M⋆ rela-
equation:
tion in Yang et al. (2018) flattens at low redshift, but ours
do not show such a trend. Therefore, the BHAR in Yang
log M⋆ = 10.65 + 0.81 log(z) + 0.83 log(1 + z), (31)
et al. (2018) is less reliable at z ≲ 0.8, as they noted; if
their relation is further extrapolated below z = 0.5, their
and they are explicitly plotted in Figure 10. We also plot
BHAR − M⋆ relation would become flat and is thus unphys-
the BHAR of the whole population in Figure 10, and it is
ical. Our BHAR uncertainties are also considerably larger
similar to the star-forming BHAR below the M⋆ threshold in
than those in Yang et al. (2018), even though we used more
Equation 31 and becomes more in the middle between the
data. This is because Yang et al. (2018) adopted a paramet-
star-forming and quiescent BHAR with rising M⋆ . There-
ric modeling method, which includes strong a priori assump-
fore, Equation 31 can also serve as an approximate threshold
tions. In contrast, this work minimizes such assumptions,
of whether the contribution of the SMBH growth in quiescent
and thus the fitted uncertainties reflect those directly inher-
galaxies to the overall SMBH growth is important.
ited from the data.
Our results suggest that the BHAR(M⋆ , z) function may
4.3. Star-forming versus Quiescent Galaxies also depend on SFR, with star-forming galaxies hosting en-
hanced AGN activity (e.g., Aird et al. 2018, 2019; Birchall
Star-forming galaxies generally have stronger AGN activ-
et al. 2023). However, such a dependence is only secondary
ity than quiescent galaxies (e.g., Aird et al. 2018, 2019). We
(Yang et al. 2017), and SFR is usually more challenging to
hence examine if star-forming and quiescent galaxies have
measure and more subject to confusion with AGN emission.
the same BHAR in this section.
Still, more physical insights can be gained by incorporating
To separate star-forming and quiescent galaxies, we adopt
SFR-based parameters, especially when probing p(λ | M⋆ , z)
the star-forming main sequence (MS) in Popesso et al. (2023)
instead of BHAR (Aird et al. 2018). We leave further analy-
and define quiescent galaxies as those lying at least 1 dex be-
ses on including SFR into the BHAR(M⋆ , z) function to the
low the MS; the remaining galaxies are star-forming ones.
future, in which different classification schemes from binary
Since the star-forming and quiescent subpopulations do not
(star-forming versus quiescent) up to four categories (star-
individually follow the SMF and XLF, we do not apply the
burst, star-forming, transitioning, and quiescent) will be ex-
SMF-XLF constraints as in Section 4.1. We measure their
plored.
BHAR and present the results in Figure 10. The BHAR of
both star-forming and quiescent galaxies increases with M⋆
and z. When comparing the BHAR of these two subpopu- 5. SUMMARY AND FUTURE WORK
lations, star-forming galaxies generally have larger BHAR, In this work, we mapped BHAR as a function of (M⋆ , z)
suggesting that star-forming galaxies indeed host more ac- over the vast majority of cosmic time, and our main results
tive SMBHs, possibly due to more available cold gas for both are summarized as follows:
18

1
z = 0.5 z = 0.8 z = 1.0
logBHAR (M yr 1)

0
1
2
3
4
1
z = 1.5 z = 2.0 z = 3.0 z = 4.0
logBHAR (M yr 1)

0
1
2
3
4
9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0
logM (M ) logM (M ) logM (M ) logM (M )
Figure 9. The comparison of our BHAR with those in Yang et al. (2018). The blue curves are our median BHAR, and the cyan shaded regions
represent the corresponding 1σ and 2σ uncertainty ranges. The BHAR and the corresponding 1σ uncertainty in Yang et al. (2018) are plotted
as the black curves.
1. We compiled an unprecedentedly large sample from those in Yang et al. (2018) at z ≳ 0.8, and we also,
nine fields – CANDELS (including GOODS-S, for the first time, provide reasonable constraints upon
GOODS-N, EGS, and UDS), the LSST DDFs (in- BHAR at lower redshift (z ≲ 0.5). See Sections 3.2
cluding COSMOS, ELAIS-S1, W-CDF-S, and XMM- and 3.3.
LSS), and eFEDS. These fields include both deep,
small-area surveys and shallow, large-area ones. The 4. We measured BHAR for both star-forming and, for
former provide rich information in the high-z, low- the first time, quiescent galaxies. Both groups show
M⋆ , and/or low-λ regime, while the latter provide com- BHAR increases with M⋆ and z, and the star-forming
plementary information in the low-z, high-M⋆ , and/or BHAR is generally larger than or at least comparable
high-λ regime. Therefore, our sample can effectively to the quiescent BHAR across the whole (M⋆ , z) plane.
constrain BHAR across a large range of the relevant See Section 4.3.
parameter space. See Section 2.
It should be noted that, besides BHAR, our p(λ | M⋆ , z)
2. We developed a semiparametric Bayesian method to parameter maps in Figure 2 also contain rich information,
measure BHAR, where a double power-law model and we release p(λ | M⋆ , z) and the corresponding parame-
with respect to λ is used to measure p(λ | M⋆ , z), and ter maps and BHAR maps in https://doi.org/10.5281/zenodo.
the relevant parameters nonparametrically depend on 10729248. As first examples, we have briefly and phe-
(M⋆ , z). This method has two main advantages. It nomenologically discussed different scientific questions in
avoids the “extrapolation” of parameters from well- Sections 4.2 and 4.3, which justified that our results can re-
populated regions in the parameter space to less- veal interesting dependences of SMBH growth on the galaxy
populated regions. It also adopts much weaker as- population.
sumptions than parametric methods, enabling more Figure 3 visually illustrates that p(λ | M⋆ , z) evolves over
flexible constraints and more representative fitting un- (M⋆ , z). Observationally, it is still unclear what the exact
certainties from the data. See Section 3.1. evolution pattern is, let alone the physics driving such an
evolution. It is also unknown from a theoretical perspective
3. We sampled p(λ | M⋆ , z) and measured BHAR at because no simulations appear to produce consistent evolu-
109.5 < M⋆ < 1012 M⊙ and z < 4. We have verified the tion patterns of p(λ | M⋆ , z) with the observed ones (e.g.,
fitting quality by comparing our model p(λ | M⋆ , z) Habouzit et al. 2022). It even complicates matters further
and the corresponding binned estimators and also by that p(λ | M⋆ , z) may evolve differently for star-forming and
comparing our inferred XLF with the observed one. quiescent galaxies, as proposed in a phenomenological sce-
We showed that BHAR increases with both M⋆ and z. nario in Aird et al. (2018). We leave detailed analyses of the
Our BHAR measurements are largely consistent with p(λ | M⋆ , z) evolution to a subsequent future work. We will
 19

1
z = 0.1 z = 0.3 z = 0.5
0
logBHAR (M yr 1)

1
2
3
4
5
1
z = 0.8 z = 1.0 z = 1.5
0
logBHAR (M yr 1)

1
2
3
4
5
1
z = 2.0 z = 3.0 z = 4.0
0
logBHAR (M yr 1)

1
2
3
4
5
9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0
logM (M ) logM (M ) logM (M )
Figure 10. BHAR for star-forming (blue) and quiescent (red) galaxies. The shaded regions represent 1σ uncertainty ranges. The black dashed
curves are the BHAR with all the galaxies included, i.e., those in Figure 6. The vertical black lines mark the maximum M⋆ values where
star-forming and quiescent galaxies can be reliably classified at the corresponding z (Cristello et al. 2024). Star-forming galaxies have larger
BHAR.
first identify the qualitative evolution pattern of the depen- ACKNOWLEDGMENTS
dence of p(λ | M⋆ , z) on M⋆ and z for different galaxy pop-
ulations and then develop a quantitative, parametric model We thank the anonymous referee for constructive sugges-
to depict the identified evolution pattern. With the clearly tions and comments. We thank Nathan Cristello, Joel Leja,
understood p(λ | M⋆ , z), we will address the following scien- and Zhenyuan Wang for helpful discussions. FZ, ZY, and
tific questions. Is the broad decline in SMBH growth below WNB acknowledge financial support from NSF grant AST-
z ≈ 1 due to the shift of accretion activity to smaller galaxies 2106990, Chandra X-ray Center grant AR1-22006X, the
or reduction of the typical λ? How large is the AGN “duty Penn State Eberly Endowment, and Penn State ACIS Instru-
cycle”, which is an integration of p(λ | M⋆ , z), in different ment Team Contract SV4-74018 (issued by the Chandra X-
galaxy populations? Does M⋆ modulate the duty cycle or ray Center, which is operated by the Smithsonian Astrophys-
modulate the typical outburst luminosity in the AGN phase? ical Observatory for and on behalf of NASA under contract
Is there any difference in the SMBH feeding in star-forming NAS8-03060). GY acknowledges funding from the Nether-
and quiescent galaxies? lands Research School for Astronomy (NOVA). The Chandra
ACIS Team Guaranteed Time Observations (GTO) utilized
were selected by the ACIS Instrument Principal Investigator,
Gordon P. Garmire, currently of the Huntingdon Institute for
X-ray Astronomy, LLC, which is under contract to the Smith-
sonian Astrophysical Observatory via Contract SV2-82024.
20

APPENDIX

γ1
A. GRADIENT OF THE POSTERIOR  λc , λ < λc
∂ ln p 

=

(A6)
∂λc  γ2 , λ > λc

This Appendix presents the gradient of our posterior in 
λc
Equation 20. We found that, at least in our case, analytical
differentiation enables a much higher computational speed 
and/or less memory compared with other differentiation algo- − ln λλ , λ < λc
 
∂ ln p 

= (A7)
 c
rithms. We thus adopt the analytically-derived gradient and ∂γ1 0, λ > λc


directly present the derivation results below.
First, the partial derivatives of I(γ, λ1 , λ2 , A, λc ; M⋆ , z) in
Equation 12 are 
0, λ < λc ,
∂ ln p 

∂I = (A8)

I
= ∂γ2 − ln λ , λ > λc
 
(A1) 
∂A A

λc

∇ ln L corresponding to Equation 16 can then be expressed

∂I γ
= I (A2) as follows.
∂λc λc
∂ ln L gal ∂I
"
= −ni j (γ1,i j , λmin , λc,i j , Ai j , λc,i j ; M⋆,i , z j )
∂Ai j ∂A
!−γ
∂I λ1 1 λ1
" ! #
A
=− ln + [erf(x1 ) + 1]
∂γ λc γ λc ∂I
#
2γ ln 10
+ (γ2,i j , λc,i j , +∞, Ai j , λc,i j ; M⋆,i , z j )
∂A
!−γ
λ2 1 λ2
" ! #
A
+ ln + [erf(x2 ) + 1]
2γ ln 10 λc γ λc nAGN
X ∂ ln p
ij
−γ
 10a− γ4bln210  "
 + (λ s ; Ai j , λc,i j , γ1,i j , γ2,i j ) (A9)
a
γ(ln 10)2 1 ∂A
! #
+
A   ln 10 − + × s=1
2γ ln 10 λc M⋆ η  λc M⋆ η γ

2b2
γ ln 10 γ ln 10
" ! !#
∂ ln L gal ∂I
"
erf x1 + − erf x2 + = −ni j (γ1,i j , λmin , λc,i j , Ai j , λc,i j ; M⋆,i , z j )
2b 2b ∂λc,i j ∂λ2
 γ ln 10 −γ
A  10a− 4b2  ∂I
− √  × + (γ1,i j , λmin , λc,i j , Ai j , λc,i j ; M⋆,i , z j )
2 πγb  λc M⋆ η  ∂λc

∂I
+ (γ2,i j , λc,i j , +∞, Ai j , λc,i j ; M⋆,i , z j )
  !2   !2 
γ ln 10 γ ln 10
exp − x1 +  − exp − x2 + ∂λ1
    

2b 2b
∂I
#
(A3) + (γ2,i j , λc,i j , +∞, Ai j , λc,i j ; M⋆,i , z j )
∂λc
nAGN
∂I X ∂ ln p
ij

= (2k − 3)× + (λ s ; Ai j , λc,i j , γ1,i j , γ2,i j ) (A10)

∂λk s=1
∂λc
!−γ
λk
(
A
[erf(xk ) + 1]
2λk ln 10 λc
∂ ln L gal ∂I
Ab λk
!−γ = −ni j (γk,i j , λmin , λc,i j , Ai j , λc,i j ; M⋆,i , z j )
∂γk,i j ∂γ
 
2
− √ exp −x
π(ln 10)2 γλk λc k
nAGN
−γ X ∂ ln p
ij

 10a− γ4bln210  + (λ s ; Ai j , λc,i j , γ1,i j , γ2,i j )

  !2  )
Ab γ ln 10 
+ √  exp − xk + ∂γk

 
π(ln 10)2 γλk  λc M⋆ η 

s=1
  
2b

(k = 1, 2) (A11)
(k = 1, 2) (A4)
The partial derivatives of ln πcont in Equation 19 are
Defining ln p(λ; A, λc , γ1 , γ2 ) as ln p(λ | M⋆ , z), its partial
derivatives are ∂ ln πcont N M (Xi−1, j + Xi+1, j − 2Xi j ) Nz (Xi, j−1 + Xi, j+1 − 2Xi j )
= + ,
∂ ln p 1 ∂Xi j σ2X σ2X
= (A5)
∂A A (A12)
 21

in which X denotes each one of (log A, log λc , γ1 , γ2 ), and we 0, λc < M LX



∂ψDP


 2
define X0 j ≡ X1 j , XNM +1, j ≡ XNM , j , Xi0 ≡ Xi1 , and Xi,Nz +1 ≡  ∂ψ
=

∂γ (γ1 , λc , M2 , A, λc ), M2 < λc <
LX LX LX (B20)
Xi,Nz to incorporate X’s at the boundary. ∂γ1 

 M1
 ∂ψ (γ1 , M1 , M2 , A, λc ), λc > LX


The gradient of the log-posterior in Equation 20 is thus

∂γ M1
X
∇ ln P = ∇ ln L + ∇ ln πcont . (A13)
field 
∂ψ
∂γ (γ2 , M1 , M2 , A, λc ), λc < M2


 LX
When transforming the parameter space, the gradient of the ∂ψDP



 ∂ψ
=

∂γ (γ2 , M1 , λc , A, λc ), M2 < λc <
LX LX LX (B21)
corresponding Jacobian should also be added.
∂γ2 



M1

B. GRADIENT OF THE POSTERIOR WITH THE

0, λc > LX


M1
SMF-XLF CONSTRAINTS ADDED
Based on Equation 26, we have
This Appendix presents the gradient of our posterior after
adding the SMF-XLF constraints in Equation 29. First, the ∂ϕL,mdl
partial derivatives of ψ(γ, M1 , M2 , A, λc ; LX ) in Equation 24 (LX , z) = δ j jz ×
∂Xi j
are
∂ψDP
∂ψ ψ (Ai jz , λc,i jz , γ1,i jz , γ2,i jz , MLB,i , MLB,i+1 ; LX ),
= (B14) ∂X
∂A A (B22)

∂ψ γψ where δ j jz = 0 (1) if j , jz ( j = jz ), and X denotes each

= (B15)
∂λc λc one of (A, λc , γ1 , γ2 ). The partial derivatives of ln LSMF−XLF
!−γ in Equation 27 are
∂ψ ∂ΓGI
" !
LX M1 M2
=A ϕ1 α1 + γ + 1, , ∂ϕL,mdl
∂ ln LSMF−XLF X XLF
!
∂γ Mc λc ∂ζ Mc Mc 1 1
= nk − (LX,k , zk ).
LX
!
M1 M2
!# ∂Xi j k
ϕL,mdl,k ϕL,obs,k ∂Xi j
− ln ΓGI α1 + γ + 1, ,
Mc λc Mc Mc (B23)
!−γ "
∂ΓGI
!
LX M1 M2
+A ϕ2 α2 + γ + 1, , The gradient of the posterior in Equation 29 is
Mc λc ∂ζ Mc Mc
! !# X
− ln
LX
ΓGI α2 + γ + 1,
M1 M2
, (B16) ∇ ln P = ∇ ln L + ∇ ln LSMF−XLF + ∇ ln πcont , (B24)
Mc λc Mc Mc field

!−γ where ∇ ln L and ∇ ln πcont were presented in Appendix A.

∂ψ A LX Mk
= (2k − 3) e− Mc ×
∂Mk Mc Mc λc
" !α1 +γ !α +γ # C. RESULTS WITHOUT EFEDS
Mk Mk 2
ϕ1 + ϕ2 (k = 1, 2) (B17) eFEDS is primarily observed through soft X-rays below
Mc Mc
2 keV, which are more prone to obscuration compared to our
where ∂Γ∂ζGI (ζ, x1 , x2 ) is the partial derivative relative to the other fields. To examine if our results are biased by this ef-
first argument of ΓGI (ζ, x1 , x2 ). The partial derivatives of fect, we try excluding eFEDS in this appendix, and the corre-
ψDP (A, λc , γ1 , γ2 , M1 , M2 ; LX ) in Equation 25 are sponding BHAR results are shown in Figure 11. There is not
any material systematic difference in the median BHAR af-
∂ψDP ψDP
= (B18) ter excluding eFEDS, and the uncertainty becomes larger in
∂A A
certain parameter ranges; e.g., the difference in width of the
shaded regions in Figure 11 is apparent at M⋆ ≈ 1010.8 M⊙
∂ψ

∂λc (γ2 , M1 , M2 , A, λc ), λc < M2
LX
and z = 0.5. The uncertainties generally grow by no more





∂ψ ∂ψ

than 60%. Therefore, no strong systematic biases are intro-
∂λc (γ2 , M1 , λc , A, λc ) + ∂λc (γ1 , λc , M2 , A, λc )



 LX LX

∂ψDP

duced by eFEDS, and eFEDS also helps constrain BHAR.


 LX ∂ψ
=
 h
− λ2 ∂M2 (γ2 , M1 , LλXc , A, λc )
∂λc 

 c This verifies that the absorption effects have been appropri-
∂ψ
 i L
+ (γ1 , LλXc , M2 , A, λc ) , M < λc < M
LX ately considered, as detailed in Section 3.1.1. Besides, given

 X
∂M


 1 2 1
that the LSST DDFs already cover 12.6 deg2 with sensitive

 ∂ψ (γ , M , M , A, λ ), λ > LX



∂λc 1 1 2 c c M1 HB data, eFEDS provides useful constraints but is not fully
(B19) dominant.
22

1
z = 0.1 z = 0.3 z = 0.5
0
logBHAR (M yr 1)

1
2
3
4
5
1
z = 0.8 z = 1.0 z = 1.5
0
logBHAR (M yr 1)

1
2
3
4
5
1
z = 2.0 z = 3.0 z = 4.0
0
logBHAR (M yr 1)

1
2
3
4
5
9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0
logM (M ) logM (M ) logM (M )
Figure 11. Comparison between BHAR with eFEDS included (red) and excluded (blue) in the fitting. The shaded regions represent 1σ
uncertainty ranges. The red curves are similar to the blue ones, and the red uncertainties are smaller than the blue ones in certain regimes,
indicating that eFEDS does not cause systematic biases and helps constrain BHAR.

D. IMPACT OF AGN-DOMINATED SOURCES ily focus on the AGN-dominated sources in the LSST DDFs
It is generally more challenging to reliably measure M⋆ and eFEDS.
from the galaxy component for sources with SEDs domi- M⋆ is largely constrained by the rest-frame near-infrared
nated by the AGN component. We assess if the less-reliable (NIR) data because the old-star emission peaks in the NIR.
M⋆ measurements for such sources have strong impact on For the purpose of assessing the M⋆ measurements, we de-
our BHAR results. It has been shown that the CANDELS fine a source to be “AGN-dominated” if its AGN component
fields are largely free from this potential issue (Aird et al. contributes > 50% of the rest-frame 1 µm light, as mea-
2018; Yang et al. 2018) due to their small solid angles, superb sured from its decomposed SED. A similar definition was
multi-wavelength coverage, and deep X-ray surveys. For the also adopted in Aird et al. (2018). About 10 − 15% of our
LSST DDFs and eFEDS, their X-ray surveys are wider and AGNs are classified as AGN-dominated. Note that this defi-
shallower, and thus a larger fraction of the detected AGNs nition significantly overlaps but is not the same as the broad-
are luminous and may dominate the SEDs. We thus primar- line AGN definition. In a general sense, broad-line AGNs
are sources with strong AGN signatures (e.g., spectroscopi-
 23

cally detected broad emission lines) in the optical. However, cluded after removing AGN-dominated sources. These pro-
a large fraction of broad-line AGNs are not necessarily AGN- cedures are conducted for the whole population as well as
dominated in the NIR because the galaxy emission usually star-forming and quiescent galaxies. We compare BHAR
reaches a peak while the AGN emission reaches a valley in with the original ones in Figure 12. The quiescent curves al-
the NIR. We found that around half of the broad-line AGNs in most do not change after removing AGN-dominated sources,
Ni et al. (2021b) are classified as AGN-dominated under our while the whole-population and star-forming BHAR become
definition, and the non-AGN-dominated ones indeed gener- slightly smaller. The difference at high redshift is slightly
ally have lower LX . We adopt our current definition because larger than that at low redshift because high-z sources need
it is simpler and also more physically related to the M⋆ mea- higher LX to be detected in the X-ray and are hence more
surement. likely to be AGN-dominated, but the difference is still gener-
We remove AGN-dominated sources in the LSST DDFs ally no more than the 1σ uncertainties. Besides, our number-
and eFEDS and measure BHAR again following Section 3.3. based correction underestimates the real loss of accretion
We further estimate the AGN number-density maps in the power because AGN-dominated sources, by construction,
(M⋆ , z) plane using kernel density estimations before and af- tend to have higher λ than the remaining ones. The differ-
ter excluding these AGN-dominated sources and apply the ence in BHAR should be even smaller. Therefore, the rela-
number-density ratio as a function of (M⋆ , z) as a correction tively larger M⋆ uncertainties of AGN-dominated sources are
of BHAR to account for the fact that fewer AGNs are in- not expected to cause material biases to our BHAR.

REFERENCES
Aird, J., Coil, A. L., & Georgakakis, A. 2017, MNRAS, 465, 3390, Brunner, H., Liu, T., Lamer, G., et al. 2022, A&A, 661, A1,
doi: 10.1093/mnras/stw2932 doi: 10.1051/0004-6361/202141266
—. 2018, MNRAS, 474, 1225, doi: 10.1093/mnras/stx2700 Chen, C. T. J., Brandt, W. N., Luo, B., et al. 2018, MNRAS, 478,
—. 2019, MNRAS, 484, 4360, doi: 10.1093/mnras/stz125 2132, doi: 10.1093/mnras/sty1036
Aird, J., Coil, A. L., & Kocevski, D. D. 2022, MNRAS, 515, 4860, Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62,
doi: 10.1093/mnras/stac2103 doi: 10.3847/0004-637X/819/1/62
Aird, J., Coil, A. L., Moustakas, J., et al. 2012, ApJ, 746, 90, Cristello, N., Zou, F., Brandt, W. N., et al. 2024, ApJ, 962, 156,
doi: 10.1088/0004-637X/746/1/90 doi: 10.3847/1538-4357/ad2177
Ananna, T. T., Treister, E., Urry, C. M., et al. 2019, ApJ, 871, 240, Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485,
doi: 10.3847/1538-4357/aafb77 4817, doi: 10.1093/mnras/stz712
Barro, G., Pérez-González, P. G., Cava, A., et al. 2019, ApJS, 243, Driver, S. P., Bellstedt, S., Robotham, A. S. G., et al. 2022,
22, doi: 10.3847/1538-4365/ab23f2
MNRAS, 513, 439, doi: 10.1093/mnras/stac472
Bayer, A. E., Seljak, U., & Modi, C. 2023, arXiv e-prints,
Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73,
arXiv:2307.09504, doi: 10.48550/arXiv.2307.09504
doi: 10.1051/0004-6361/201936817
Betancourt, M. 2017, arXiv e-prints, arXiv:1701.02434,
Feldmann, R. 2017, MNRAS, 470, L59,
doi: 10.48550/arXiv.1701.02434
doi: 10.1093/mnrasl/slx073
Birchall, K. L., Watson, M. G., & Aird, J. 2020, MNRAS, 492,
Gehrels, N. 1986, ApJ, 303, 336, doi: 10.1086/164079
2268, doi: 10.1093/mnras/staa040
Georgakakis, A., Aird, J., Schulze, A., et al. 2017, MNRAS, 471,
Birchall, K. L., Watson, M. G., Aird, J., & Starling, R. L. C. 2023,
1976, doi: 10.1093/mnras/stx1602
MNRAS, 523, 4756, doi: 10.1093/mnras/stad1723
Bongiorno, A., Merloni, A., Brusa, M., et al. 2012, MNRAS, 427, Georgakakis, A., Nandra, K., Laird, E. S., Aird, J., & Trichas, M.
3103, doi: 10.1111/j.1365-2966.2012.22089.x 2008, MNRAS, 388, 1205,
Bongiorno, A., Schulze, A., Merloni, A., et al. 2016, A&A, 588, doi: 10.1111/j.1365-2966.2008.13423.x
A78, doi: 10.1051/0004-6361/201527436 Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS,
Brandt, W. N., & Alexander, D. M. 2015, A&A Rv, 23, 1, 197, 35, doi: 10.1088/0067-0049/197/2/35
doi: 10.1007/s00159-014-0081-z Habouzit, M., Somerville, R. S., Li, Y., et al. 2022, MNRAS, 509,
Brandt, W. N., & Yang, G. 2022, in Handbook of X-ray and 3015, doi: 10.1093/mnras/stab3147
Gamma-ray Astrophysics, 78, Hickox, R. C., Mullaney, J. R., Alexander, D. M., et al. 2014, ApJ,
doi: 10.1007/978-981-16-4544-0 130-1 782, 9, doi: 10.1088/0004-637X/782/1/9
Brandt, W. N., Ni, Q., Yang, G., et al. 2018, arXiv e-prints, Jarvis, M. J., Bonfield, D. G., Bruce, V. A., et al. 2013, MNRAS,
arXiv:1811.06542, doi: 10.48550/arXiv.1811.06542 428, 1281, doi: 10.1093/mnras/sts118
24

1
z = 0.1 Included z = 0.3 z = 0.5
0 Excluded
logBHAR (M yr 1)

1
2
3
4
5
1
z = 0.8 z = 1.0 z = 1.5
0
logBHAR (M yr 1)

1
2
3
4
5
1
z = 2.0 z = 3.0 z = 4.0
0
logBHAR (M yr 1)

1
2
3
4
5
9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0 9.5 10.0 10.5 11.0 11.5 12.0
logM (M ) logM (M ) logM (M )
Figure 12. Comparison between BHAR with AGN-dominated sources included (solid curves) and excluded (dashed curves) in the fitting.
Black, blue, and red curves represent the whole population, star-forming galaxies, and quiescent galaxies, respectively. The grey shaded
regions are 1σ uncertainty ranges of the black solid curves. The solid and dashed BHAR curves are generally consistent within 1σ uncertainties,
indicating that AGN-dominated sources do not cause material biases to our results.

Kocevski, D. D., Hasinger, G., Brightman, M., et al. 2018, ApJS, Loredo, T. J. 2004, in American Institute of Physics Conference
236, 48, doi: 10.3847/1538-4365/aab9b4 Series, Vol. 735, Bayesian Inference and Maximum Entropy
Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, Methods in Science and Engineering: 24th International
ApJS, 197, 36, doi: 10.1088/0067-0049/197/2/36 Workshop on Bayesian Inference and Maximum Entropy
Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, Methods in Science and Engineering, ed. R. Fischer, R. Preuss,
doi: 10.1146/annurev-astro-082708-101811 & U. V. Toussaint, 195–206, doi: 10.1063/1.1835214
Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24, Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2,
doi: 10.3847/0067-0049/224/2/24 doi: 10.3847/1538-4365/228/1/2
Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, Lusso, E., Comastri, A., Simmons, B. D., et al. 2012, MNRAS,
J. S. 2019, ApJ, 876, 3, doi: 10.3847/1538-4357/ab133c 425, 623, doi: 10.1111/j.1365-2966.2012.21513.x
Liu, T., Buchner, J., Nandra, K., et al. 2022, A&A, 661, A5, Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34,
doi: 10.1051/0004-6361/202141643 doi: 10.3847/0004-637X/817/1/34
 25

Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, ApJS, 220, 10, Xue, Y. Q., Luo, B., Brandt, W. N., et al. 2016, ApJS, 224, 15,
doi: 10.1088/0067-0049/220/1/10 doi: 10.3847/0067-0049/224/2/15
Ni, Q., Yang, G., Brandt, W. N., et al. 2019, MNRAS, 490, 1135, Yan, W., Brandt, W. N., Zou, F., et al. 2023, ApJ, 951, 27,
doi: 10.1093/mnras/stz2623 doi: 10.3847/1538-4357/accea6
Ni, Q., Brandt, W. N., Yang, G., et al. 2021a, MNRAS, 500, 4989,
Yang, G., Brandt, W. N., Alexander, D. M., et al. 2019, MNRAS,
doi: 10.1093/mnras/staa3514
485, 3721, doi: 10.1093/mnras/stz611
Ni, Q., Brandt, W. N., Chen, C.-T., et al. 2021b, ApJS, 256, 21,
doi: 10.3847/1538-4365/ac0dc6 Yang, G., Estrada-Carpenter, V., Papovich, C., et al. 2021, ApJ,
Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 921, 170, doi: 10.3847/1538-4357/ac2233
1526, doi: 10.1093/mnras/stac3214 Yang, G., Brandt, W. N., Luo, B., et al. 2016, ApJ, 831, 145,
Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, doi: 10.3847/0004-637X/831/2/145
A13, doi: 10.1051/0004-6361/200913020 Yang, G., Chen, C. T. J., Vito, F., et al. 2017, ApJ, 842, 72,
Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes
doi: 10.3847/1538-4357/aa7564
for Machine Learning
Yang, G., Brandt, W. N., Vito, F., et al. 2018, MNRAS, 475, 1887,
Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017, Nature, 549,
doi: 10.1093/mnras/stx2805
488, doi: 10.1038/nature23906
Salvato, M., Wolf, J., Dwelly, T., et al. 2022, A&A, 661, A3, Yang, G., Caputi, K. I., Papovich, C., et al. 2023, ApJL, 950, L5,
doi: 10.1051/0004-6361/202141631 doi: 10.3847/2041-8213/acd639
Santini, P., Ferguson, H. C., Fontana, A., et al. 2015, ApJ, 801, 97, Yu, Z., Zou, F., & Brandt, W. N. 2023, Research Notes of the
doi: 10.1088/0004-637X/801/2/97 American Astronomical Society, 7, 248,
Stefanon, M., Yan, H., Mobasher, B., et al. 2017, ApJS, 229, 32, doi: 10.3847/2515-5172/ad0ed7
doi: 10.3847/1538-4365/aa66cb
Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529,
Tak, H., Ghosh, S. K., & Ellis, J. A. 2018, MNRAS, 481, 277,
doi: 10.1146/annurev-astro-082812-141003
doi: 10.1093/mnras/sty2326
Zou, F., Yang, G., Brandt, W. N., et al. 2021, Research Notes of the
Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G.
2014, ApJ, 786, 104, doi: 10.1088/0004-637X/786/2/104 American Astronomical Society, 5, 56,
Wang, T., Elbaz, D., Alexander, D. M., et al. 2017, A&A, 601, doi: 10.3847/2515-5172/abf050
A63, doi: 10.1051/0004-6361/201526645 Zou, F., Brandt, W. N., Chen, C.-T., et al. 2022, ApJS, 262, 15,
Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, doi: 10.3847/1538-4365/ac7bdf
11, doi: 10.3847/1538-4365/ac3078 Zou, F., Brandt, W. N., Ni, Q., et al. 2023, ApJ, 950, 136,
Wright, A. H., Driver, S. P., & Robotham, A. S. G. 2018, MNRAS, doi: 10.3847/1538-4357/acce39
480, 3491, doi: 10.1093/mnras/sty2136