# TRINITY - III. Quasar luminosity functions decomposed by halo, galaxy, and black hole masses as well as Eddington ratios from z = 0-10

Haowen Zhang, Peter Behroozi, Marta Volonteri, Joseph Silk, Xiaohui Fan, James Aird, Jinyi Yang, Philip F. Hopkins

► To cite this version:

Haowen Zhang, Peter Behroozi, Marta Volonteri, Joseph Silk, Xiaohui Fan, et al. TRINITY - III. Quasar luminosity functions decomposed by halo, galaxy, and black hole masses as well as Eddington ratios from z = 0-10. Monthly Notices of the Royal Astronomical Society, 2024, 529, pp.2777-2793. 10.1093/mnras/stae655. insu-04851569

HAL Id: insu-04851569

https://insu.hal.science/insu-04851569v1

Submitted on 20 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# TRINITY – III. Quasar luminosity functions decomposed by halo, galaxy, and black hole masses as well as Eddington ratios from $z = 0–10$

Haowen Zhang (张昊文)<fnref n="1" />,<fnref n="1" />★† Peter Behroozi<fnref n="1" />,<fnref n="1" />,<fnref n="2" /> Marta Volonteri<fnref n="3" /> Joseph Silk⁹,<fnref n="3" />,<fnref n="4" />,<fnref n="5" /> Xiaohui Fan<fnref n="1" />
James Aird⁹,<fnref n="6" />,<fnref n="7" /> Jinyi Yang (杨锦怡)<fnref n="1" /> and Philip F. Hopkins⁹

Accepted 2024 March 1. Received 2024 January 22; in original form 2023 June 1

## ABSTRACT

We present the redshift evolution of quasar luminosity functions (QLFs) decomposed by halo mass, galaxy mass, supermassive black hole (SMBH) mass, and Eddington ratio, as well as SMBH/radiative energy output ratios from TRINITY, a flexible empirical model that self-consistently infers the halo–galaxy–SMBH connection that match observational data. Key findings include: (1) The normalization of quasar luminosity function increases by ~3–4 dex from $z \sim 10$ to $z \sim 4$, due to the fast mass build-up of different SMBH populations; (2) From $z \sim 4$ to $z \sim 1$, less massive galaxies and SMBHs make up bigger and bigger fractions of QLFs, due to the active galactic nuclei (AGNs) downsizing effect; (3) At $z \sim 0$, massive haloes/galaxies/SMBHs are responsible for most bright quasars due to low Eddington ratios among all SMBHs; (4) The bright ends of QLFs are dominated by SMBHs that are at least 0.3 dex overmassive relative to the median SMBH mass–galaxy mass relation; (5) QLFs at $z \sim 6–7$ are dominated by SMBHs accreting at Eddington ratios $0.1 < \eta_{rad} < 1$, but super-Eddington AGNs contribute more significantly to QLFs towards $z \sim 9–10$.

**Key words:** galaxies: evolution – galaxies: haloes – quasars: supermassive black holes.

## 1 INTRODUCTION

Supermassive black holes (SMBHs) are found in the centres of many galaxies (Kormendy & Richstone 1995; Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Ho 2008; Gültekin et al. 2009; Kormendy & Ho 2013; Heckman & Best 2014). When actively accreting, SMBHs release huge amounts of energy in radiative and/or kinetic forms. This energy can regulate SMBHs' own growth and the star formation of their host galaxies (e.g. Silk & Rees 1998; Bower et al. 2006; Somerville et al. 2008; Sijacki et al. 2015). On the other hand, galactic physical processes such as supernova feedback and galaxy mergers can also regulate SMBH growth by, e.g. modulating gas fueling and/or leading to SMBH mergers. Hence, galaxy–SMBH interactions (and potential coevolution) are a critical component of galaxy formation and evolution theories (see e.g. Hopkins et al. 2007; Ho 2008; Alexander & Hickox 2012; Kormendy & Ho 2013; Heckman & Best 2014; Brandt & Alexander 2015).

Different lines of evidence support a coevolution scenario. First, there are tight galaxy–SMBH scaling relations (e.g. the black hole mass–velocity dispersion relation or black hole mass–bulge mass relation) in the local Universe (see Häring & Rix 2004; Gültekin et al. 2009; Kormendy & Ho 2013; McConnell & Ma 2013; Savorgnan et al. 2016). Second, studies have found a roughly constant ratio of $10^{-4}–10^{-3}$ between the cosmic SMBH accretion rate (CBHAR) density and the cosmic star formation rate (CSFR) density (Merloni, Rudnick & Di Matteo 2004; Silverman et al. 2008; Shankar, Weinberg & Miralda-Escudé 2009; Aird et al. 2010; Delvecchio et al. 2014; Yang et al. 2018; Peca et al. 2023; Zhang et al. 2023c) over $0 < z < 4$. However, it has been traditionally hard to examine other predictions of coevolution, like a tight SMBH–galaxy mass connection at higher redshifts. This is due to: (1) systematic uncertainties in SMBH and galaxy mass measurements at high redshifts (see Shen 2013 and references therein, also Izumi et al. 2019); (2) difficulty in measuring masses for non-active SMBHs, galaxy mass, and velocity dispersion at $z > 0$, due to the limits in sensitivity and resolution (see the summary by, e.g. Reines & Volonteri 2015); and (3) the bias that overmassive SMBHs are overrepresented in luminosity-selected samples (Lauer et al. 2007).

<footnote n="1">¹Department of Astronomy, University of Arizona, 933 N Cherry Avenue, Tucson, AZ 85721, USA</footnote>

<footnote n="2">²Division of Science, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="3">³Institut d'Astrophysique de Paris (UMR 7095: CNRS & Sorbonne Universite), 98 bis Bd. Arago, Paris F-75014, France</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA</footnote>

<footnote n="5">⁵BIPAC, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="6">⁶Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="7">⁷Department of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK</footnote>

<footnote n="8">⁸Theoretical Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote>* E-mail: hwzhang0595@arizona.edu</footnote>

<footnote>† Strittmatter Fellow.</footnote>
Several theoretical models have invoked kinetic energy as the primary driver of SMBH–galaxy coevolution (e.g. Croton et al. 2006; Somerville et al. 2008; Fanidakis et al. 2012; Hirschmann et al. 2012; Costa, Sijacki & Haehnelt 2014; Volonteri et al. 2016; Choi et al. 2017; Weinberger et al. 2017). When SMBHs are actively accreting, gravitational energy from accreted matter can be converted into both radiation and kinetic energy. At high Eddington ratios, geometrically thin but optically thick accretion discs form around AGNs, which converts accreted matter into radiation effectively (see e.g. Shakura & Sunyaev 1973). At low Eddington ratios, accretion flows becomes radiatively inefficient, and converts more accreted matter into kinetic energy (see theoretical studies like Narayan & Yi 1994 and observational studies like Nagar, Falcke & Wilson 2005, Ho 2002 and references therein). Hence, theoretical models often assumed that the fraction of kinetic energy output dominates in the low Eddington ratio regime (Merloni & Heinz 2008). However, the measurement of AGN kinetic energy is often indirect, relying on observations of X-ray cavities in gaseous haloes (Hlavacek-Larrondo et al. 2015; McDonald et al. 2021) or loose scaling relations between kinetic energy and radio luminosities (Merloni & Heinz 2007; Cavagnolo et al. 2010; Smolčić et al. 2017). As a result, the kinetic energy fraction is less well quantified when SMBHs are accreting at higher Eddington ratios. Shankar et al. (2008) assumed a simple linear relation between the AGN kinetic and radiative luminosities, with scatter. By fitting to the AGN optical and radio luminosity functions from Cirasuolo, Magliocchetti & Celotti (2005), they obtained a median kinetic-to-radiative luminosity ratio of $0.1^{+0.05}_{-0.01}$, and a log-normal scatter of $0.38^{+0.04}_{-0.09}$. Based on the scaling relations between radio luminosities (the proxy for AGN kinetic energy output), X-ray luminosities (the proxy for radiative energy output), and black hole masses, Merloni & Heinz (2008) estimated AGN kinetic versus radiative energy outputs as functions of redshift. La Franca, Melini & Fiore (2010) took this a step further by modelling the full probability distribution of AGN radio luminosity as a function of X-ray luminosity and redshift, and obtained the cosmic AGN kinetic power density as a function of redshift.

This paper builds on many past studies that have estimated the buildup of the SMBH population and quasar luminosity functions (QLFs). Beyond the local Universe, the buildup of the SMBH population is studied mostly by observing actively accreting SMBHs, i.e. active galactic nuclei (AGNs). The total cosmic mass in SMBHs can grow only via accretion, even though individual SMBHs can grow via mergers. Based on the Soltan argument (Soltan 1982), the cosmic SMBH mass density ($\rho_{BH}(z)$) is the integral of the quasar luminosity function ($\Phi(L, z)$) along the dimensions of luminosity and cosmic time, modulated by the radiative efficiency, $\epsilon_{rad}$ of AGNs:

$$ \rho_{BH}(z) = \int_{0}^{z} \int_{0}^{\infty} \frac{1 - \epsilon_{rad}}{\epsilon_{rad} c^2} L \Phi(L, z) dL dt . \quad (1) $$

To calculate the temporal and luminosity-weighted average $\epsilon_{rad}$, one compares the integral with the local SMBH black hole mass density $\rho_{BH}(z = 0)$. This is traditionally found by convolving the distribution function of a certain galaxy property $X_{gal}$ and the $X_{gal}$-SMBH mass scaling relation, with the implicit assumption that SMBHs exist in the centres of all galaxies. Based on different AGN luminosity functions, different studies generally find average $\epsilon_{rad} \sim 0.1$ (Yu & Tremaine 2002; Marconi et al. 2004; Merloni et al. 2004; Shankar et al. 2009; Shankar, Weinberg & Miralda-Escudé 2013).

In TRINITY, we apply the Soltan argument to different halo/galaxy populations, rather than treating all the galaxies in the Universe as a single population. This enables us to decompose QLFs by many halo/galaxy/SMBH properties, and helps us better understand: (1)

The contributions from different halo/galaxy/SMBH populations to total QLFs; (2) Typical AGN activity levels throughout the host galaxy growth histories; (3) How selection effects impact observed galaxy-SMBH mass connections at $z > 0$; and (4) The amplitude of AGN clustering. We also predict QLFs beyond $z \sim 6$, and investigate the dominant specific SMBH accretion rates and plausible formation mechanisms of SMBHs with $M_{\bullet} \gtrsim 10^9 M_{\odot}$ at high redshifts ($z \gtrsim 6$).

In this paper, we continue to present findings from TRINITY, an empirical model that infers self-consistent halo–galaxy–SMBH connections from $z = 0 - 10$. Compared to previous empirical models, TRINITY is equipped with a larger compilation of galaxy and AGN data. To constrain the halo–galaxy connection, we included galaxy stellar mass functions (SMFs), galaxy UV luminosity functions (UVLFs), galaxy quenched fractions (QFs), galaxy specific star formation rates (SSFRs), and cosmic star formation rates (CSFRs). The information on the galaxy–SMBH connection comes from QLFs, quasar probability distribution functions (QPDFs), active black hole mass functions (ABHMFs), and the local bulge mass–black hole mass relations. The inclusion of QPDFs is crucial, because they constrain AGN duty cycles and SMBH growth histories of different halo/galaxy/SMBH populations. Several previous studies made assumptions about Eddington ratio distributions, such as a mass-independent Eddington ratio that may change with redshift (Marconi et al. 2004; Shankar et al. 2009; Caplar, Lilly & Trakhtenbrot 2015, 2018), or a fixed Eddington ratio distribution (Shankar et al. 2013). With these assumptions, AGN duty cycles and SMBH growth trajectories can be obtained. None the less, these models represent a simplified picture of the Universe, as Eddington ratios are observed to vary with both redshift and host stellar mass (e.g. Aird, Coil & Georgakakis 2018).

The joint constraining power of this unique data set allows TRINITY to have a more flexible parametrization than past studies, which includes, e.g. the redshift evolution and host mass dependence of AGN Eddington distribution shapes, and duty cycles.

TRINITY inherits the empirical modelling framework of Behroozi, Wechsler & Conroy (2013). TRINITY makes a guess for the halo–galaxy–SMBH connection over time and constructs a mock universe by applying this guess to simulated halo populations. This mock universe is compared with the real Universe in terms of the observables above, and a Markov Chain Monte Carlo (MCMC) algorithm is used to obtain the posterior distribution of allowed halo–galaxy–SMBH connections. Based on this redshift-dependent halo–galaxy–SMBH connection, Zhang et al. (2023a) (Paper I hereafter) reconstructed average SMBH growth histories as a function of host halo/galaxy mass, and found that: (1) the median $M_{\bullet}-M_{*}$ relation experiences minimal evolution for galaxies above $5 \times 10^9 M_{\odot}$, and (2) SMBH activity experiences downsizing, in the sense that more massive SMBHs experience decline in Eddington ratios earlier than less massive SMBHs. Zhang et al. (2023b) (Paper II hereafter) utilized the modelled random scatter around the $M_{\bullet}-M_{*}$ relation and SMBH Eddington ratio distributions to quantify the luminosity-dependent bias in the observed $M_{\bullet}-M_{*}$ relation for $z \sim 6$ AGNs (i.e. Lauer bias, Lauer et al. 2007).

As the third in the TRINITY paper series, this work discusses QLFs and the buildup of supermassive black holes across cosmic time. The fourth paper (Paper IV) makes predictions for SMBHs at $z > 6$. The fifth paper (Paper V) investigates the variation of the SFR–BHAR correlation with halo/galaxy/SMBH mass and redshift. In the sixth paper (Paper VI) we calculate black hole merger statistics and make predictions for gravitational wave instruments including NANOGray and LISA. The seventh (Paper VII) and eighth (Paper VIII) papers present the AGN autocorrelation functions and
AGN-galaxy cross-correlation functions from TRINITY, respectively. They also discuss whether/how models like TRINITY can be constrained by AGN clustering signals.

The paper is organized as follows. Section 2 covers methodology. In Section 3, we describe the dark matter simulations and galaxy/SMBH observations used to constrain TRINITY. Section 4 presents the results of our model, followed by the discussion and the comparison with other models in Section 5. Finally, we present the conclusions in Section 6. In this work, we adopt a flat $\Lambda$CDM cosmology with parameters ($\Omega_m = 0.307$, $\Omega_\Lambda = 0.693$, $h = 0.678$, $\sigma_8 = 0.823$, $n_s = 0.96$) consistent with Planck results (Planck Collaboration XIII 2016). We use data sets that adopt the Chabrier stellar initial mass function (IMF; Chabrier 2003), the Bruzual & Charlot (2003) stellar population synthesis model, and the Calzetti dust attenuation law (Calzetti et al. 2000). Halo masses are calculated following the virial overdensity definition from Bryan & Norman (1998).

## 2 METHODOLOGY

In Section 2.1 we explain: (1) why TRINITY is able to infer the halo-galaxy-SMBH connection from $z = 0–10$ uniquely, robustly, and self-consistently; and (2) why such inference can be done by fitting observations alone, without assuming any particular physical mechanisms for, e.g. halo-galaxy-SMBH interactions. Section 2.2 gives an overview of the TRINITY implementation, and the calculation of AGN observable properties are laid out in Section 2.3.

### 2.1 Why observations alone can constrain the halo-galaxy-SMBH connection

The Soltan argument (Soltan 1982) has traditionally been applied as a cosmic average: the ratio of the total luminosity output of SMBHs to their $z=0$ mass density gives the cosmic average radiative efficiency, and this in turn allows inferring the cosmic average growth history of SMBHs from the redshift evolution of the total luminosity in QLFs.

More recently, studies including Shankar et al. (2020) have recognized that the Soltan argument applies for any sub-population of SMBHs. That is, the average growth history of SMBHs in a chosen set of galaxies at $z=0$ is proportional to the net AGN luminosity of those galaxies’ progenitors at earlier times. Shankar et al. (2020) chose to estimate SMBH histories in bins of host galaxy stellar mass by combining AGN luminosity measurements as a function of galaxy mass and redshift with average galaxy growth histories (i.e. approximating each galaxy as having only a single progenitor at earlier times). In TRINITY, we forward model AGN luminosities across the distribution of galaxy growth histories to: (1) improve the treatment of scatter in galaxy accretion histories and mergers with other galaxies; (2) self-consistently include selection effects; and (3) self-consistently incorporate additional galaxy and SMBH data sources.

In both Shankar et al. (2020) and this paper, a well-constrained halo-galaxy relation is essential to inferring the halo-galaxy-SMBH connection, as well as galaxy assembly histories. Fortunately, such tight constraints are provided in recent empirical models of the halo-galaxy connection, e.g. Moster, Naab & White (2013, 2018) and Behroozi et al. (2013); Behroozi et al. (2019). In these models, the galaxy-halo relation is constrained by the number density and clustering of galaxies as a function of redshift. Since halo formation histories are well-quantified in dark matter simulations, a constraint on the halo-galaxy connection automatically gives a constraint on galaxy star formation histories. Incorrect galaxy formation histories

(e.g. formation times for Milky Way-mass galaxies that are too short) lead to tension with observations (e.g. number densities for Milky Way progenitor galaxies that are too low compared to observations, as well as specific star formation rates that are too high compared with observations) – meaning that constraints on galaxy formation histories are strong and data-driven.

Last but not least, the robustness of TRINITY predictions is well tested. We experimented with an extensive series of model variants by changing: (1) the parametrization of the galaxy-SMBH mass connection; (2) accretion rate distribution shapes; (3) SMBH merger prescriptions; (4) AGN obscuration corrections; (5) AGN bolometric corrections, and found that our predictions do not change qualitatively when the input observations are self-consistent. For full details, we refer readers to the Appendices of Zhang et al. (2023c).

### 2.2 Implementation overview

TRINITY is an empirical model that infers the statistical connection between dark matter haloes, galaxies, and supermassive black holes from $z=0–10$ (Zhang et al. 2023c). Unlike the UNIVERSEMA-CHINE, which follows individual haloes across cosmic time, TRINITY keeps track of typical galaxy and SMBH masses of different halo populations as functions of redshift. In this section, we summarize the main aspects of the model. For full details about the TRINITY parametrization, we refer readers to section 2 of Zhang et al. (2023c). For clarity, ‘host haloes’ of galaxies and SMBHs in this work (and in all TRINITY papers) refers to individual subhaloes that host single galaxies at their centres, rather than the parent haloes that host multiple satellite haloes/galaxies.

TRINITY connects dark matter haloes and galaxies by parametrizing the average galaxy SFR as a function of peak historical host halo mass ($M_{peak}$)<fnref n="1" /> and redshift. Any guess for the SFR($M_{peak}$, $z$) function allows populating a simulated halo population with galaxy SFRs across cosmic time. By integrating these SFRs over halo assembly histories, we establish the stellar mass-halo mass ($M_* - M_{peak}$) relation as a function of redshift. At a given halo mass and redshift, galaxy stellar mass is converted into bulge mass based on a redshift-dependent bulge mass-total stellar mass ($M_{bulge} - M_*$) relation from SDSS (Mendel et al. 2014) and CANDELS (Lang et al. 2014) observations. With galaxy bulge masses in place, median and average black hole masses are assigned to haloes with a parametrized black hole mass-bulge mass ($M_\bullet - M_{bulge}$) relation that is allowed to evolve with time. Every combination of the $M_*$-$M_{peak}$ and $M_\bullet$-$M_{bulge}$ relations thus fully specifies the average galaxy stellar and SMBH masses given the host halo mass and redshift. Although SMBH occupation fraction is formally a free parameter in TRINITY, we find in practice that all haloes with $M_{peak} > 10^{11} M_\odot$ are constrained to host SMBHs to match $z=0$ constraints (see also Habouzit, Volonteri & Dubois 2017; Tremmel 2017; Zhang et al. 2023c). Specifically, QLFs and quasar probability distribution functions (QPDFs) constrain the total amount of SMBH mass growth for different galaxy mass bins. A small occupation fraction will lead to too few SMBHs to account for this mass growth, resulting in a small number of very massive BHs rather than many less massive BHs. As such, a small occupation fraction will overproduce low-redshift active black hole mass functions (ABHMFs) at the massive end.

Average total SMBH growth rates at fixed halo mass are calculated by taking the difference of SMBH masses assigned to halo

<footnote n="1">¹In detail, TRINITY parametrizes SFR($v_{Mpeak}(M_{peak}, z)$, $z$), where $v_{Mpeak}$ is the maximum circular velocity when the halo reaches its peak mass, $M_{peak}$.</footnote>
populations between two consecutive snapshots. Given continuing uncertainties in SMBH merger rates, we assume that the fractional merger contribution to total SMBH growth is proportional to the fractional merger contribution to total galaxy growth, with a redshift-dependent proportionality factor.

To convert average SMBH accretion rates to observable AGN luminosity distributions, we parametrize the AGN total energy efficiency ($\epsilon_{tot}$, including radiative+kinetic energies, redshift-independent), duty cycles, and Eddington ratio distribution shapes. Many of our observational constraints are on AGN radiative power, so we adopt a scaling relation between the radiative ($\epsilon_{rad}$), kinetic ($\epsilon_{kin}$), and total ($\epsilon_{tot} \equiv \epsilon_{rad} + \epsilon_{kin}$) energy efficiencies based on previous studies (Mineshige et al. 2000; Merloni & Heinz 2008). The scaling between $\epsilon_{rad}$ and $\epsilon_{kin}$ is a function of SMBH Eddington ratio:

$$ \eta = \frac{\epsilon_{tot} \text{BHAR}/(M_{\odot}/\text{yr}) \times 4.5 \times 10^8}{(1 - \epsilon_{tot}) M_{\odot}/M_{\odot}}, \quad (2) $$

where BHAR is SMBH accretion rate. This scaling relation is inspired by the observational study by Merloni & Heinz (2007). We refer the readers to Appendix A for the full details of the $\epsilon_{rad} - \epsilon_{kin}$ scaling relation.

To account for the systematic uncertainties induced by our assumed $\epsilon_{rad} - \epsilon_{tot}$ relation, we also reran TRINITY with two different alternative models: (1) the one where no kinetic energy is produced from SMBH accretion (the ‘no kinetic’ model hereafter), and (2) the one where more kinetic energy is generated in SMBH accretion (the ‘more kinetic’ model hereafter). These two variants represent the two extremes of kinetic energy output in SMBH accretion, but our main conclusions remain invariant when these alternative models are adopted. Full details are introduced in Appendix A. We assume that SMBHs in the same halo mass bin share the same Eddington ratio distribution.

Given that observational studies usually calculate Eddington ratios based on radiative luminosities rather than total energy output, we define ‘kinetic’ and ‘radiative’ Eddington ratios:

$$ \eta_{\text{kin/rad}} = \frac{\epsilon_{\text{kin/rad}}(\eta)}{\epsilon_{\text{tot}}} \eta. \quad (3) $$

Systematic effects are also included as nuisance parameters in this work, including systematic offsets in stellar mass and SFR measurements with SED fitting, and the systematic offsets in Eddington ratios when predicting different AGN observables. We also adopt empirical corrections from observational studies for other systematic effects such as Compton-thick AGN obscuration fractions. After correcting for systematic effects, we forward model the joint distribution of galaxy and AGN observables and compare it with real observations (summarized in section 3 of Zhang et al. 2023c). Such a comparison yields the likelihood of the model parameters given the real data constraints. To calculate the posterior distribution of the model parameters, we adopt a custom Metropolis Markov Chain Monte Carlo (MCMC) algorithm (based on Haario, Saksman & Tamminen 2001). The Markov chain length is set to 2 million steps to ensure convergence, which is approximately 50 autocorrelation lengths. From this posterior distribution, we obtain both: (1) the best-fitting halo–galaxy–SMBH connection; and (2) the uncertainties in this connection allowed by current observations. The best-fitting model parameters, down-sampled MCMC chain, and correlations between model parameters can be found at https://github.com/HaowenZhang/TRINITY/tree/main/stats.

## 2.3 Calculating black hole statistics

In this work, we discuss the build-up of SMBHs. This build-up process is reflected in different statistics, whose calculations are covered here.

TRINITY effectively parameterizes AGN bolometric luminosity ($L_{bol}$), and black hole mass ($M_{\bullet}$) distributions as functions of halo mass and redshift (i.e. $P(L_{bol}|M_{\bullet}, M_{peak}, z)$ and $P(M_{\bullet}|M_{peak}, z)$), from which all other observables are computed. The cosmic density of a certain AGN quantity, X, at each redshift is calculated by adding up the contributions from all the haloes in the Universe:

$$ \rho_X(z) = \int_0^\infty \bar{X}(M_{\text{peak}}, z) \phi_h(M_{\text{peak}}, z) dM_{\text{peak}}, \quad (4) $$

where ‘X’ can be either black hole accretion rate (BHAR) or bolometric luminosity ($L_{bol}$, Section 4.3), etc., and $\phi_h(M_{peak}, z)$ is the dark matter halo mass function at redshift z.

QLFs (Section 4.2) are given by the number density of haloes hosting SMBHs with a given luminosity:

$$ \phi_L(L_{\text{bol}}, z) = \int_0^\infty \phi_h(M_{\text{peak}}, z) P(L_{\text{bol}}|M_{\text{peak}}, z) dM_{\text{peak}}, \quad (5) $$

where $P(L_{bol}|M_h, z)$ is calculated by counting the number density of black holes with the corresponding Eddington ratio:

$$ P(L_{\text{bol}}|M_{\text{peak}}, z) = \int_0^\infty P(L_{\text{bol}}, |M_{\bullet}, M_{\text{peak}}, z) \times P(M_{\bullet}|M_{\text{peak}}, z) dM_{\bullet}, \quad (6) $$

The AGN duty cycle (Section 4.4) can be defined in several different ways. If defined as the fraction of black holes above a certain luminosity threshold, $L_{min}$, then the duty cycle $f(> L_{min})$ is:

$$ f(> L_{\min} | M_{\text{peak}}, z) = \int_{L_{\min}}^{\infty} P(L_{\text{bol}} | M_{\text{peak}}, z) dL_{\text{bol}}. \quad (7) $$

Similarly, if defined as the fraction of black holes above a certain Eddington ratio threshold, $\eta_{min}$, then the duty cycle $f(> \eta_{min})$ is:

$$ \begin{aligned} f(> \eta_{\min} | M_{\text{peak}}, z) &= \int_{\eta_{\min}}^{\infty} P(\eta_{\text{rad}} | M_{\text{peak}}, z) d\eta_{\text{rad}} \\ &= \int_{\eta_{\min}}^{\infty} P(\eta | M_{\text{peak}}, z) \frac{d\eta}{d\eta_{\text{rad}}} d\eta_{\text{rad}} \end{aligned} \quad (8) $$

where $P(\eta|M_{peak}, z)$ is the (conditional) Eddington ratio distribution as a function of host halo mass and redshift.

# 3 SIMULATIONS AND DATA CONSTRAINTS

## 3.1 Dark matter halo statistics

Unlike Behroozi et al. (2019), TRINITY does not follow individual haloes over cosmic time, but only uses halo statistics from N-body simulations. Specifically, we use halo mass functions from Behroozi et al. (2013) for the cosmology specified in the introduction. These mass functions are similar to the central halo mass functions from Tinker et al. (2008), but small adjustments were made to account for satellite haloes as well as to use $M_{peak}$ (the historical peak halo mass) rather than the current mass. These adjustments were based on the Bolshoi & Consuelo simulations (Klypin, Trujillo-Gomez & Primack 2011). The full details of these halo mass functions can be found in appendix G of Behroozi et al. (2013). Note that the Bolshoi simulation was based on the WMAP7 cosmology (Jarosik et al. 2011) rather than the Planck cosmology, which is used everywhere else in TRINITY. As such, we made further corrections to these halo mass
<table><caption>**Table 1.** Summary of observational constraints.</caption><thead><tr><th>Type</th><th>Redshifts</th><th>Primarily constrains</th></tr></thead><tbody><tr><td>Stellar mass functions</td><td>0–8</td><td>SFR–M<sub>peak</sub> relation</td></tr><tr><td>Galaxy quenched fractions</td><td>0–4</td><td>Quenching–M<sub>peak</sub> relation</td></tr><tr><td>Cosmic star formation rates</td><td>0–10</td><td>SFR–M<sub>peak</sub> relation</td></tr><tr><td>Specific star formation rates</td><td>0–9</td><td>SFR–M<sub>peak</sub> relation</td></tr><tr><td>Galaxy UV luminosity functions</td><td>9–10</td><td>SFR–M<sub>peak</sub> relation</td></tr><tr><td>QLFs</td><td>0–5</td><td>Total SMBH accretion</td></tr><tr><td>Quasar probability distribution functions</td><td>0–2.5</td><td>AGN duty cycle, BHAR distributions</td></tr><tr><td>Active black hole mass functions</td><td>0–5</td><td>Scatter around the galaxy–SMBH mass connection</td></tr><tr><td>Black hole mass—bulge mass relation</td><td>0</td><td>Galaxy–SMBH connection, AGN energy efficiency</td></tr><tr><td>Observed black hole mass distribution of bright quasars</td><td>5.8–6.5</td><td>Galaxy–SMBH connection, Eddington ratios at high redshifts</td></tr></tbody></table>

*Note.* BHAR is the black hole accretion rate.

functions to account for the difference between the cosmologies. This correction factor is calculated by comparing the total (i.e. central+satellite) halo mass functions from Behroozi et al. (2013) to those from the Bolshoi-Planck simulation (Rodríguez-Puebla et al. 2016), as a function of halo mass and redshift. With this correction, these calibrated halo mass functions are suitable for studying the evolution of haloes from $10^{10}$ to $10^{15}$ M$_{\odot}$. Of note, the halo mass functions used here are for individual subhaloes that host single galaxies in their centres, rather than the parent haloes that host multiple galaxies and/or SMBHs.

Both accretion and mergers contribute to halo growth. We use the fitting formulae in appendix H of Behroozi et al. (2013) to describe the average halo accretion histories. As for halo mergers, the merger rates are fitted from the UNIVERSEMACHINE (Behroozi et al. 2019). The fitting formulae for halo mergers are available in appendix B of Zhang et al. (2023c).

## 3.2 Observational data constraints

To constrain the halo–galaxy connection, we have compiled the following galaxy observables: stellar mass functions (SMFs), galaxy quenched fractions (QFs), average specific star formation rates (SSFRs), cosmic star formation rates (CSFRs), and galaxy UV luminosity functions (UVLFs). Collectively, these data cover a redshift range of $z = 0–10$. The galaxy–SMBH connection is constrained by the following AGN observables: QLFs (from Ueda et al. 2014), quasar probability distribution functions (QPDFs, from Aird et al. 2018), active black hole mass functions (ABHMFs, from Schulze & Wisotzki 2010; Kelly & Shen 2013; Schulze et al. 2015), the local $M_{\bullet}$–$M_{\text{bulge}}$ relation, and the observed $M_{\bullet}$ distribution of high redshift bright quasars. The redshift coverage of these AGN data is from $z = 0–6.5$. All the data sets are summarized in Table 1, and we refer readers to Behroozi et al. (2019) and Zhang et al. (2023c) for the full details about galaxy and AGN observables, respectively.

JWST has detected more high-redshift UV-bright galaxies compared to the predictions of many theoretical models, e.g. the UNIVERSEMACHINE. This discrepancy between new observations and existing model predictions means that these latest data provide new information to constrain the high-redshift behaviours of our theoretical models. In light of this, we added the $9 \le z \le 13$ spectroscopically confirmed galaxy UVLFs from Harikane et al. (2024) to our data compilation. In TRINITY, we calculate galaxy UV magnitudes with a scaling relation between UV magnitudes and SFR, to get statistically equivalent results to the Flexible Stellar Population Synthesis code (FSPS; Conroy & White 2013). The full details of calculating galaxy UV magnitudes can be found in appendix B of Zhang et al. (2023c).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The comparison of luminosity functions between the original TRINITY model (dashed lines, Zhang et al. 2023c), the updated TRINITY model (solid lines), the existing galaxy UVLF data (stars, Ishigaki et al. 2018; Oesch et al. 2018; Bouwens et al. 2019, also see table 8 of Zhang et al. 2023c), and the new galaxy UVLFs from JWST (Harikane et al. 2024). The difference between the original and updated TRINITY models increases with redshift. See Section 3.2.</figcaption></figure>

Adding $9 \le z \le 13$ galaxy UVLFs does not induce any qualitative changes in the results presented by Zhang et al. (2023a, c). Quantitatively, this addition leads to more UV-bright galaxies at $z \gtrsim 9$ and thus increased galaxy mass at a given halo mass. The redshift-dependent $M_{\bullet}$–$M_{\star}$ relation is unchanged with newly added galaxy UVLFs, and thus there will be more SMBH mass at fixed halo mass as well. Fig. 1 shows the change in galaxy number density, which increases towards higher redshifts. This is because the new JWST UVLFs (solid circles) are largely consistent with existing $z \le 10$ UVLFs in our data compilation (solid stars) But the original TRINITY underproduced UV-bright galaxies without the higher redshift UVLFs, leading to the larger difference at these redshifts. The results presented in this work are based on the updated TRINITY model after adding these $9 \le z \le 13$ galaxy UVLFs to our data constraints.

We note that since the SMBH data only cover a redshift range of $z = 0–6.5$, all the predictions for $z \gtrsim 6$ in this work are effectively extrapolations based on lower redshift data constraints and our parametrization. Future observations by, e.g. JWST, Roman, or Euclid will provide critical tests of our $z \gtrsim 6$ predictions. Recently, JWST has detected many new SMBHs at $z \gtrsim 6$. We have compared some of the JWST detections in another paper focused on $z \gtrsim 6$ predictions (Paper IV), and will defer more thorough comparisons to future studies when the systematic and selection effects in observations are better quantified.
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** The evolution of the bolometric quasar luminosity function from 0.1 ≤ z ≤ 2.5 (top panel) and from 2.5 ≤ z ≤ 10 (bottom panel). The solid parts of the curves represent regions where QLF data is available, and the dashed parts denote the predictions beyond the redshift or luminosity ranges covered by observed QLFs. The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. See Section 4.1.</figcaption></figure>

# 4 RESULTS

Section 4.1 shows the redshift evolution of QLFs from z = 0–10. Unless otherwise noted, all the QLFs shown in this work include both non-obscured and obscured (including Compton-thick obscured) quasars. In Section 4.2, we present the QLFs split into different bins of SMBH mass, galaxy mass, halo mass, deviation from the intrinsic $M_{\bullet}-M_{*}$ relation, and Eddington ratio. AGN radiative power densities across cosmic time are shown in Section 4.3. AGN duty cycles from TRINITY are presented in Section 4.4. We then show the fraction of super-Eddington AGNs in Section 4.5. For the comparisons of the best-fitting model to observations, we refer readers to Zhang et al. (2023c).

## 4.1 Redshift evolution of QLFs from z = 0–10

The top panel of Fig. 2 shows the evolution of quasar bolometric luminosity function (QLF) from 0.1 ≤ z ≤ 2.5. Towards z ~ 2.5, there is increasing AGN activity in the universe. We thus see the decrease(increase) in the faint(bright) end of QLF, and an increase in the typical quasar luminosity from z = 0.1 to z = 2.5. These trends are qualitatively consistent with Shen et al. (2020). The bottom panel shows bolometric QLF at 2.5 ≤ z ≤ 10. In Fig. 2, the solid parts of the curves represent regions where QLF data is available, whereas the dashed parts denote the predictions beyond the redshift

or luminosity ranges covered by observed QLFs. As was found by Shen et al. (2020), the normalization of QLF drops quickly towards higher redshifts at z > 4, by ~3–4 dex from z ~ 4 to z ~ 10. This is simply due to the lack of massive SMBHs to drive bright AGNs in the early universe. We also see increasing uncertainties around the QLF at higher redshifts, because black hole accretion rates at higher redshifts are smaller, have less impact on the model predictions at lower redshifts, and thus are not as well constrained by existing observations.

## 4.2 Quasar luminosity function decompositions

### 4.2.1 Quasars binned by $M_{\bullet}$

Fig. 3 shows quasar bolometric luminosity functions (QLFs) in bins of black hole mass at different redshifts. There are kinks at the faint end of QLFs for different mass bins. They are due to the change in the scaling of $\eta_{\text{tot}}$ with $\eta_{\text{tot}}$ at $\eta_{\text{tot}} = 0.03$ (see Appendix A).

At lower redshifts, e.g. z ≲ 1, the QLFs of different black hole populations decrease monotonically with luminosity; at higher redshifts, however, the QLFs start to show a peak. This results from the narrowing of Eddington ratio distributions towards higher redshifts. These higher and narrower Eddington ratio distributions are required by the combination of QPDFs from Aird et al. (2018) and the ABHMFs from Kelly & Shen (2013). Specifically, the QPDFs specifies the average SMBH accretion histories in different galaxies (by specifying the average SMBH accretion rates via numerical integrals) as well as AGN duty cycles (i.e. fraction of SMBHs accreting at certain luminosities). With such constraints, lower Eddington ratios and/or broader Eddington ratio distributions cannot reproduce the ABHMFs at high-z. At 6 ≲ z ≲ 8 (9 ≲ z ≲ 10), all black holes in haloes above ~$10^{11.5} M_{\odot}$ are accreting at around (above) the Eddington rate. As a result, different black hole mass populations dominate the QLF at different luminosities.

### 4.2.2 Quasars binned by host galaxy and host halo mass

Fig. 4 shows QLFs in bins of galaxy stellar mass at different redshifts. At z ≳ 6, galaxies of different masses dominate different parts of the QLFs. This is again due to the fact that all active SMBHs in haloes above ~$10^{11.5} M_{\odot}$ are accreting at around the Eddington rate (see Section 4.2.1). The spreads in luminosity for each stellar mass bin are broader than those for each bin in SMBH mass. This is due to: 1) the higher-than-unity $M_{\bullet}-M_{*}$ slope at high redshifts; and (2) scatter in the $M_{\bullet}-M_{*}$ relation. From z ~ 4 to z ~ 1, galaxies with $10 < \log_{10}M_{*} < 11$ become more and more prevalent compared to those with $11 < \log_{10}M_{*} < 12$. This is mainly due to the stronger decline in the activity level of SMBHs in more massive galaxies with time, i.e. the AGN downsizing effect. From z ~ 1 to z ~ 0.1, galaxies with $10 < \log_{10}M_{*} < 11$ lose dominance in the bright end of QLFs. At z = 0.1, galaxies with $9 < \log_{10}M_{*} < 11$ dominate QLFs below $L_{\text{bol}} \sim 3 \times 10^{45} \text{ erg s}^{-1}$; above this luminosity threshold, more massive objects with $11 < \log_{10}M_{*} < 12$ are dominant. This over-representation of massive galaxies at the bright end of the QLF is primarily due to the similar or higher Eddington ratios for high-mass SMBHs/haloes than their lower-mass counterparts, as required by the low-redshift QPDFs from Aird et al. (2018). These $M_{*}$ predictions for different quasar populations will help us better understand the expected host galaxy properties of different quasars in current and future observations, e.g. with JWST and ALMA.

Fig. 5 shows the average host galaxy masses ($\langle M_{*} \rangle$) for AGNs above different bolometric luminosity limits as functions of redshift.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** QLFs in bins of black hole mass at $z = 0.1$ (upper left), $z = 1.0$ (upper right), $z = 3.0$ (middle left), $z = 6.0$ (middle right), $z = 8.0$ (lower left), and $z = 10.0$ (lower right). The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. See Section 4.2.1.</figcaption></figure>

In general, $\langle M_* \rangle$ increases with cosmic time from $z = 10$ to $z \sim 5$, where active SMBHs are growing at high Eddington ratios with their host galaxies. $\langle M_* \rangle$ then peak and decrease as the ‘AGN downsizing’ effect becomes significant: SMBHs in more massive haloes/galaxies become less active, and those in smaller galaxies take over to drive AGNs. Below $z \lesssim 0.5$, the AGN activity level further decreases, such that bright AGNs must be driven by massive enough SMBHs. As a result, we see the increase in $\langle M_* \rangle$ towards $z = 0$, which is the more significant for brighter AGN populations. We also show the 68th percentile range of host $M_*$ for AGNs with $\log L_{bol} \ge 43$

(lighter shaded region) and $\log L_{bol} \ge 47$ (darker shaded region). Since only massive galaxies can host bright AGNs, the $M_*$ range for bright AGNs is in general narrower than that for fainter AGNs.

Fig. 6 shows QLFs in bins of host halo<fnref n="2" /> mass at different redshifts. At $z \gtrsim 6$, haloes with $10^{11} < M_{\text{peak}} < 10^{12} M_{\odot}$ ($10^{12} < M_{\text{peak}} < 10^{13} M_{\odot}$) dominate the faint (bright) end of QLFs. The critical

<footnote n="2">²The ‘host halo’ here refers to individual subhaloes that host single galaxies in their centres, rather than the parent haloes that host multiple galaxies and/or SMBHs.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** QLFs in bins of galaxy stellar mass at $z = 0.1$ (first row, left), $z = 0.2$ (first row, right), $z = 0.3$ (second row, left), $z = 1.0$ (second row, right), $z = 4.0$ (third row, left), $z = 6.0$ (third row, right), $z = 8.0$ (bottom row, left), and $z = 10.0$ (bottom row, right). The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. See Section 4.2.2.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The average host galaxy masses ($\langle M_* \rangle$) for AGNs above different bolometric luminosity limits, as functions of redshift. The lighter(darker) shaded region is the 68th percentile range of host $M_*$ for AGNs with log $L_{bol} \ge 43(47)$. See Section 4.2.2.</figcaption></figure>

luminosity where the contributions from these two halo mass bins are comparable decreases towards lower redshifts. This is because at higher redshifts, there are fewer massive haloes to dominate QLFs, except for the brightest end. From $z=3$ to $z=1$, haloes with $10^{12} < M_{\text{peak}} < 10^{13} M_{\odot}$ dominate the whole QLF. Towards lower redshifts, they are gradually taken over by lower mass haloes at the faint end, largely due to the AGN downsizing effect. From $z=1$ to $z=0.1$, $10^{12} < M_{\text{peak}} < 10^{13} M_{\odot}$ haloes are joined by $M_{\text{peak}} > 10^{13} M_{\odot}$ haloes at the bright end. This is constrained by the local ABHMF from Schulze & Wisotzki (2010) and the local $M_{\bullet}-M_{\text{bulge}}$ relation, which require that massive AGNs live in $M_{\text{peak}} > 10^{13} M_{\odot}$ haloes.

The relative importance of $M_{\text{peak}}\sim 10^{12} M_{\odot}$ haloes in the QLF arises because these haloes are most efficient at forming galaxies, which in turn results in the most black hole growth compared to haloes in other mass ranges. As a result, studies that probe quasar clustering usually find similar inferred host halo masses, with relatively weak redshift and luminosity trends. This conclusion is further discussed in Papers VII and VIII, which focus on AGN clustering.

### 4.2.3 Quasars binned by offset from the $M_{\bullet}-M_{*}$ relation

Fig. 7 shows QLFs in bins of offset in SMBH mass compared to the median $M_{\bullet}-M_{*}$ relation at different redshifts. This figure quantifies the Lauer bias, i.e. the amount by which brighter quasars are more likely to be driven by over-massive SMBHs (compared to their host galaxy mass) than typical SMBHs (Lauer et al. 2007). Observations constrain this effect mainly because extremely massive host galaxies are very rare by themselves, and cannot account for the observed number density of bright quasars (e.g. Yang et al. 2021, also see Zhang et al. 2023a). According to TRINITY, overmassive SMBHs ($>0.3$ dex above the median $M_{\bullet}-M_{\text{bulge}}$ relation) contribute significantly to QLFs down to $\log_{10}L_{\text{bol}} \sim 46–47$ beyond $z \sim 1$. At $z \sim 0$, this limit is lower, $\log_{10}L_{\text{bol}} \sim 45.3$. This results from the overall decline in SMBH activity level. Typical SMBHs around the intrinsic $M_{\bullet}-M_{*}$ relation dominate QLFs at $\log_{10}L_{\text{bol}} \sim 45$. Such a strong Lauer bias demonstrates the necessity of including faint AGN (i.e. $\log_{10}L_{\text{bol}} \sim 45$) when studying the redshift evolution of the galaxy–SMBH mass connection, instead of comparing the local $M_{\bullet}-M_{*}$ relation with quasar mass estimates. Undermassive SMBHs are always subdominant contributors to QLFs throughout cosmic time.

Given that Zhang et al. (2023b, c) have already shown the $M_{\bullet}-M_{*}$ relations for luminosity-selected AGNs at $z \ge 6$, we refer readers to these two papers for the plots of such scaling relations at $z \ge 6$, and https://github.com/HaowenZhang/TRINITY/tree/main/plots/Paper3/Additional_BHSM_lbol for the plots at $z < 6$.

### 4.2.4 Quasars binned by eddington ratio

Fig. 8 shows QLFs in Eddington ratio bins at different redshifts. At low redshifts, Eddington ratio distributions are broader. Consequently, quasars with different Eddington ratios dominate different parts of QLFs, with higher Eddington ratio objects being generally brighter. Between $1 \le z \le 3$, super-Eddington AGNs dominate the QLFs at $\log L_{\text{bol}} \gtrsim 47.5$. This comes from the extrapolation of QPDFs from Aird et al. (2018) towards highest specific black hole accretion rates (sBHARs). In TRINITY, we chose not to include QPDFs from Aird et al. (2018) at the high-sBHAR end. This is because they adopted a smoothness prior on modeled QPDFs, which may lead to systematic issues when the sample size is small. A larger sample of massive galaxies with highly active SMBHs is needed for more robust determination of QPDFs in this regime, which will serve as a strong test of this prediction by TRINITY. We also note that the MCMC uncertainties for QLFs in different Eddington ratio bins are small between $1 \le z \le 3$. This is because the parametrization of the TRINITY model is not flexible enough to cover the full uncertainty ranges as determined by Aird et al. (2018). Consequently, our estimates of QLF uncertainties in different Eddington ratio bins are likely lower limits. Towards higher redshifts, SMBHs of all masses become more active. In other words, Eddington ratio distributions become narrower and shift towards higher values. At $6 \le z \le 8$, essentially all bright SMBHs are accreting at close to the Eddington rate. As a result, high-redshift QLFs are dominated by quasars with $0.1 < \eta_{\text{rad}} < 1.0$. At $9 \le z \le 10$, the contribution from super-Eddington quasars becomes more significant. The recent detection and measurement of the SMBH in GN-z11 is consistent with this picture (Maiolino et al. 2023), although a larger sample from future observations by, e.g. JWST, is needed for a statistical comparison.

## 4.3 AGN radiative power densities

The top panel of Fig. 9 shows the cosmic AGN radiative power densities, $\rho_{\text{rad}}$, from TRINITY and Merloni & Heinz (2008). Qualitatively, $\rho_{\text{rad}}$ (including the total $\rho_{\text{rad}}$ from different SMBH populations) from TRINITY is broadly consistent with that from Merloni & Heinz (2008) below $z \sim 2$. But black holes in TRINITY have much higher Eddington ratios at $z \gtrsim 4$. As a result, we also predict higher radiative power densities for every mass bin compared to Merloni & Heinz (2008) in the early Universe. Another significant difference is in the low-redshift behaviour of black holes with $10^6 M_{\odot} < M_{\bullet} < 10^7 M_{\odot}$: Merloni & Heinz (2008) predicts a power density that increases with time below $z \sim 0.5$, while we find the opposite. One potential explanation is that in addition to luminosity functions, we included the quasar probability distribution functions from Aird et al. (2018). These data put upper limits on the energy output from low mass black holes in TRINITY, which were absent in Merloni & Heinz (2008). The bottom panel of Fig. 9 shows the same $\rho_{\text{rad}}$, as well as the contributions from different galaxy populations. In terms of redshift evolution, $\rho_{\text{rad}}$ from different galaxy populations are qualitatively similar to those of different SMBH populations, including the ‘AGN downsizing’ effect. There are offsets between curves of the same colour from these two panels, especially for lower
<figure><img src="image_9.png" /><figcaption>**Figure 6.** QLFs in bins of host halo mass at $z = 0.1$ (upper left), $z = 1.0$ (upper right), $z = 3.0$ (middle left), $z = 6.0$ (middle right), $z = 8.0$ (lower left), $z = 10.0$ (lower right). The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. See Section 4.2.2.</figcaption></figure>

mass galaxies/SMBHs at higher redshifts. This is mostly driven by the stronger redshift evolution in the $M_\bullet-M_*$ relation at the less-massive end.

## 4.4 Black hole duty cycles

TRINITY predicts SMBH Eddington ratio distributions as functions of halo/galaxy/SMBH mass and redshift. This enables us to calculate AGN duty cycles based on observables such as Eddington ratios or

luminosities. Fig. 10 shows AGN duty cycles when defined as the fractions of black holes accreting at $\eta_{\text{rad}} > 0.01$ ($f(\eta_{\text{rad}} > 0.01)$, top panel) and $\eta_{\text{rad}} > 0.1$ ($f(\eta_{\text{rad}} > 0.1)$, bottom panel), as a function of $M_{\text{peak}}$ and $z$. In general, both fractions increase with redshift. In general, $f(\eta_{\text{rad}} > 0.01)$ and $f(\eta_{\text{rad}} > 0.1)$ are significantly lower in haloes with $\log M_{\text{peak}} \lesssim 11.5$. This is constrained by the QPDFs from Aird et al. (2018), which indicates lower duty cycles for low-mass galaxies. Note that these QPDFs only cover a redshift range of $z \le 2.5$. Thus, the prediction for low-mass haloes beyond $z \sim 2.5$
<figure><img src="image_10.png" /><figcaption>**Figure 7.** QLFs in bins of the deviation in SMBH mass from the median $M_{\bullet}-M_{*}$ relation at $z = 0.1$ (upper left), $z = 1.0$ (upper right), $z = 3.0$ (middle left), $z = 6.0$ (middle right), $z = 8.0$ (lower left), and $z = 10.0$ (lower right). The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. See Section 4.2.3.</figcaption></figure>

are effectively extrapolations by TRINITY. Overall, such low AGN duty cycles in low-mass haloes mean that SMBHs in these low-mass haloes/galaxies likely grow in intermittent accretion events. The majority of their growth histories are instead dominated by a dormant/low-activity phase. During this dormant/low-activity phase, any AGN feedback for these dwarf galaxies is likely dominated by kinetic energy, as suggested in Appendix A. These low duty cycles in low-mass haloes is also seen in many recent theoretical studies, e.g.

Dubois et al. 2015; Anglés-Alcázar et al. 2017; Bower et al. 2017; Habouzit et al. 2017 (also see the summary by Tillman et al. 2022). In these models, the gas fuel for SMBH growth are likely removed by supernova feedback in low-mass galaxies. While their host galaxies keep growing, these SMBHs cannot grow substantially via extended and intensive accretion events. SMBHs growth becomes efficient after their host galaxies are too massive for supernova feedback to be efficient, which allows gas reservoir to build up in galaxy
<figure><img src="image_11.png" /><figcaption>**Figure 8.** QLFs split in different Eddington ratio slices at $z = 0.1$ (upper left), $z = 1.0$ (upper right), $z = 3.0$ (middle left), $z = 6.0$ (middle right), $z = 8.0$ (lower left) and $z = 10.0$ (lower right). The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. The red shaded regions denotes the typical luminosity range for observed $z \sim 6$ quasars. See Section 4.2.4.</figcaption></figure>

centers. The low AGN duty cycles can also solve a common issue in many theoretical models of galaxy formation: many of them are not calibrated against QLFs. So they tend to overproduce the faint ends of QLFs (e.g. Sijacki et al. 2015; Volonteri et al. 2016; Habouzit et al. 2017; Habouzit et al. 2022). By matching both the faint-end QLFs from Ueda et al. (2014) and the QPDFs for low-mass galaxies from Aird et al. (2018), TRINITY shows that it is possible to reproduce the faint-end QLFs by keeping the majority of low-mass

SMBHs dormant across cosmic time, and only allow intermittent AGN activity.

Above $z \sim 5$, both $f(\eta_{\text{rad}} > 0.01)$ and $f(\eta_{\text{rad}} > 0.1)$ are essentially mass-independent, above $\log M_{\text{peak}} \sim 11.5$. This is because the Eddington ratio distributions of all SMBHs become narrower and shifted to higher values at higher redshifts. Below $z \sim 3$, both fractions decrease towards lower redshifts. More massive black holes experience this decline earlier, known as the ‘AGN downsizing’
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Top panel: Cosmic AGN radiative power densities from TRINITY and Merloni & Heinz (2008). Different colours denote the contribution from black holes of various masses. The shaded regions around the solid curves are the 68 per cent confidence intervals from TRINITY, while other shaded regions denote the predictions from Merloni & Heinz (2008). Bottom panel: The same as the top panel, but different colours denote the contributions from different stellar mass bins. See Section 4.3.</figcaption></figure>

phenomenon. Note that there is residual AGN activity among massive haloes with $M_{\text{peak}} \gtrsim 10^{13.5} M_{\odot}$ between $0.5 \le z \le 2$. This is extrapolated from the QPDFs from Aird et al. (2018), which implies higher average SMBH accretion rate among more massive galaxies. Overall, the decrease in Eddington ratio also coincides with the decline of the star-forming fraction among galaxies, as shown by the $f_{\text{SF}} = 0.5$ curve (yellow dashed line). Specifically, the halo mass at which the fractions of star-forming and quiescent galaxies become equal (i.e. the $f_{\text{SF}} = 0.5$ curve) coincides with an average Eddington ratio of $\eta_{\text{rad}} \sim 0.01$.

We also show AGN duty cycles when defined as the fraction of black holes with bolometric luminosities $L_{\text{bol}}$ above a certain limit in Fig. 11. The top panel shows the fraction of black holes with $L_{\text{bol}} > 10^{43} \text{ erg s}^{-1}$ ($f_{>43}$), and the bottom panel shows the fraction of black holes with $L_{\text{bol}} > 10^{45} \text{ erg s}^{-1}$ ($f_{>45}$), as functions of $M_{\text{peak}}$ and $z$. The mass and redshift dependencies of both fractions are similar to the ones in Fig. 10. But, compared to low mass black holes, there are more massive black holes with $L_{\text{bol}} > 10^{43} \text{ erg s}^{-1}$ than $\eta_{\text{rad}} > 0.01$ below $z \sim 3$. The underlying reason is simple: luminosities depend on the product of black hole mass and Eddington ratio. Thus, despite lower Eddington ratios, massive black holes can still produce similar luminosities as smaller objects. In general, AGN duty cycles

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Top panel: the fractions of black holes accreting at $\eta_{\text{rad}} > 0.01$ as a function of $M_{\text{peak}}$ and $z$. Bottom panel: the fractions of black holes accreting at $\eta_{\text{rad}} > 0.1$ as a function of $M_{\text{peak}}$ and $z$. The white solid lines are the average mass growth curves of haloes with $M_{\text{peak}} = 10^{12}$, $10^{13}$, $10^{14}$, and $10^{15} M_{\odot}$ at $z=0$. The yellow dashed line shows the halo mass at which the galaxy star-forming fraction $f_{\text{SF}}$ is 0.5 as a function of $z$. See Section 4.4.</figcaption></figure>

are higher at high redshifts. One implication of this is that the active black hole mass functions (ABHMFs) become closer to total BHMFs (i.e. including both active and dormant SMBHs) at high redshifts.

## 4.5 Super-eddington accretion

The derivation of the Eddington limit means that it represents at most a time-averaged limit on SMBH accretion for spherically symmetric systems. Indeed, one of the findings of Zhang et al. (2023c) was that existing SMBH observations require super-Eddington accretion to occur some fraction of the time. Specifically, Eddington-limited accretion yields too slow SMBH growth to reproduce the local $M_{\bullet}-M_{\text{bulge}}$ relation (see appendix E1 of Zhang et al. 2023c). On the other hand, super-Eddington accretion is also necessary to reproduce sufficient numbers of bright quasars and SMBH mass functions simultaneously.

Fig. 12 shows the fractions of black holes accreting at super-Eddington rates as a function of $M_{\text{peak}}$ and $z$. The mass and redshift dependence of the super-Eddington fraction is qualitatively similar to that shown in Fig. 10. This demonstrates that the decline in black hole activity below $z \sim 6$ happens among black holes at all different activity levels.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Top panel: the fractions of black holes with $L_{bol} > 10^{43}$ erg s$^{-1}$ as a function of $M_{peak}$ and z. Bottom panel: the fractions of black holes with $L_{bol} > 10^{45}$ erg s$^{-1}$ as a function of $M_{peak}$ and z. The white solid lines are the average mass growth curves of haloes with $M_{peak} = 10^{12}$, $10^{13}$, $10^{14}$, and $10^{15}$ M$_{⊙}$ at z = 0. See Section 4.4.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** The fractions of black holes accreting above the Eddington rate as a function of $M_{peak}$ and z. The white solid lines are the average mass growth curves of haloes with $M_{peak} = 10^{12}$, $10^{13}$, $10^{14}$, and $10^{15}$ M$_{⊙}$ at z = 0. See Section 4.5.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Comparison between AGN efficiencies ($\epsilon_{tot}$) from different studies. The black error bar is the 10–90th percentile range from Fanidakis et al. (2012). The blue dot and error bar are the best-fitting value and the 68 per cent confidence intervals for TRINITY. The red error bar shows the range of systematic change in $\epsilon_{tot}$ when relevant input assumptions are changed (for details about this systematic change, see appendix E2 of Zhang et al. 2023c). See Section 5.1.</figcaption></figure>

## 5 COMPARISON WITH PREVIOUS STUDIES AND DISCUSSION

Section 5.1 compares the mass-to-energy conversion efficiency of AGNs in TRINITY and other studies, and Section 5.2 discusses the potential correlation between the decline in AGN Eddington ratio and the transition between AGN feedback modes.

### 5.1 AGN energy efficiency

Fig. 13 shows the inferred or assumed AGN *total* radiative + kinetic efficiencies, $\epsilon_{tot}=\epsilon_{rad} + \epsilon_{kin}$, in the following studies: TRINITY, the empirical models by Marconi et al. (2004), Merloni & Heinz (2008), Shankar et al. (2013), and Caplar et al. (2018), and the semi-analytical models by Fanidakis et al. (2012) and Hirschmann et al. (2012). Fanidakis et al. (2012) calculated energy efficiencies of individual black holes from their spins, so here we show the 10–90th percentile range that covers their black holes. For TRINITY, the blue error bar denotes the 16–84th percentile confidence interval from the posterior distribution of the model parameters, and the red error bar shows the range of systematic change in $\epsilon_{tot}$ when relevant input assumptions are changed, as discussed in appendix E2 of Zhang et al. (2023c). The $\epsilon_{tot}$ from TRINITY lies within the range spanned by different studies, but it is also important to keep in mind that the uncertainties in input assumptions to TRINITY leads to additional systematic uncertainty of ~0.25 dex in $\epsilon_{tot}$. Such uncertainties include the choice of local $M_{\bullet}-M_{bulge}$ relation, the choice of bolometric corrections, the Compton-thick fraction corrections, and the way that the galaxy–SMBH mass connection is parametrized. For further details, we refer readers to appendix E2 of Zhang et al. (2023c).

### 5.2 AGN feedback mode versus eddington ratio

In the top panel of Fig. 10, we showed that the following regions are coincident with each other: (1) the region where the fraction of star-forming (quenched) galaxies is exactly 50 per cent (SF/Q equality
region hereafter); (2) the region where the fraction of SMBHs with Eddington ratio $\eta_{\text{rad}} > 0.01$ decreases dramatically (AGN boundary hereafter). Observations require the coincidence of these regions because the stellar mass at which galaxies transition between star-forming and quiescent coincides with strong decreases in specific quasar luminosity distributions from Aird et al. (2018). Interestingly, $\eta_{\text{rad}} \sim 0.01$ is also consistent with the boundary between ‘radiative mode’ and the ‘kinetic mode’ feedback – i.e. the boundary between primarily radiative and primarily kinetic energy output from AGNs. This two-mode feedback scenario is widely adopted in semi-analytical models and hydrodynamical simulations (e.g. Bower et al. 2006; Croton et al. 2006; Sijacki et al. 2007, 2015; Somerville et al. 2008; Dubois, Volonteri & Silk 2014; Vogelsberger et al. 2014; Weinberger et al. 2017; Davé et al. 2019). However, it is also worth noting that the fraction of active SMBHs defined with respect to other critical Eddington ratio values (e.g. $\eta_{\text{rad}} = 0.1$, the bottom panel of Fig. 10) also drops sharply close to the SF/Q equality region, so it remains unclear if $\eta_{\text{rad}} \sim 0.01$ is where the actual transition between radiative mode and kinetic mode feedback occurs. In addition, this overlap of boundaries only indicates *correlation*, rather than *causation*, between SMBH quiescence and galaxy quenching.

# 6 CONCLUSIONS

In this work, we present the build-up of SMBH populations from TRINITY, which is an empirical model that parametrizes the statistical halo–galaxy–SMBH connection (Section 2). Compared to previous studies that are typically focused on one or two kinds of observables, TRINITY features the ability to self-consistently match a comprehensive set of observational data for galaxies and SMBHs from $z = 0–10$. As a result, TRINITY is able to extract information that is only accessible with joint constraints from multiple datasets. In this paper, we investigated quasar demographics; key results are as follows:

(i) From $z \sim 10$ to $z \sim 4$, the normalization of quasar bolometric luminosity functions (QLFs) increases by $\sim 3–4$ dex (Section 4.1, Fig. 2). This is due to the mass buildup of different SMBH populations in TRINITY, which is informed by observed QLFs and quasar probability distribution functions (QPDFs) at $z \le 5$ and extrapolations at $z \ge 6$.

(ii) Overall, less massive galaxies and SMBHs become more and more dominant in QLFs towards lower redshifts. This is due to the AGN downsizing effect, where more massive SMBHs become less active at earlier cosmic times. An exception is that at $z \sim 0$, massive haloes/galaxies/SMBHs dominate the bright end of QLFs again. This is largely due to the low Eddington ratios among all SMBHs, which makes massive SMBHs the only population able to support bright AGNs (Section 4.2.1 and 4.2.2 and Figs 3 and 4). These results are informed by observed QPDFs and QLFs.

(iii) At $z \gtrsim 1$, most observed quasars live in dark matter haloes with $M_{\text{peak}} \sim 10^{12} M_{\odot}$ (Section 4.2.2, Fig. 6). This is because $M_{\text{peak}} \sim 10^{12} M_{\odot}$ haloes are where the bulk of the star formation happens in the Universe, as required by galaxy data constraints, and because SMBH growth closely tracks galaxy growth in these galaxies.

(iv) Many observed quasars have SMBH masses $\gtrsim 0.3$ dex higher than predicted by the median $M_{\bullet}-M_{*}$ relation for all SMBHs, leading to significant Lauer bias in inferred $M_{\bullet}-M_{*}$ relations. At brighter luminosities, QLFs are more dominated by overmassive SMBHs compared to their host galaxies. Due to the overall decline in the AGN activity level with time, these overmassive SMBHs dominate QLFs down to lower luminosities at lower redshifts (Section 4.2.3,

Fig. 7). This magnitude of Lauer bias is determined by the Eddington ratio distribution and the random scatter around the $M_{\bullet}-M_{*}$ relation, which are constrained by the shapes of QPDFs, QLFs, and active black hole mass functions (ABHMFs).

(v) At $z \sim 6–8$, QLFs are dominated by AGNs with radiative Eddington ratios in the range $0.1–1$, i.e. mildly sub-Eddington to Eddington rates. At $z \sim 9–10$, we predict that super-Eddington AGNs contribute significantly to total QLFs (Section 4.2.4, Fig. 8). This result is driven by extrapolating SMBH growth histories towards $z = 10$ and the observational constraint that most bright quasars at $z \sim 6$ are sub-Eddington.

(vi) Below $z \sim 6$, the decline of AGN activity happens at all levels, i.e. the fraction of SMBHs above all Eddington ratio thresholds decrease towards lower redshifts (Section 4.4 and 4.5 and Figs 10 and 12). This result is driven by the redshift evolution of QLFs, QPDFs, and ABHMFs.

# ACKNOWLEDGEMENTS

We thank the anonymous referee for their valuable and constructive comments, as well as Gurtina Besla, Haley Bowden, Arjun Dey, Sandy Faber, Dan Foreman-Mackey, Richard Green, Jenny Greene, Melanie Habouzit, Andrew Hearin, Raphael Hviding, Takuma Izumi, David Koo, Pavel Kroupa, Tod Lauer, Junyao Li, Jianwei Lyu, Alessandro Peca, Joel Primack, Yujing Qin, George Rieke, Marcia Rieke, Jan-Torge Schindler, John Silverman, Xuejian Shen, Yue Shen, Rachel Somerville, Fengwu Sun, Wei-Leong Tee, Yoshihiro Ueda, Feige Wang, and Minghao Yue for very valuable discussions.

Support for this research came partially via program number HST-AR-15631.001-A, provided through a grant from the Space Telescope Science Institute under NASA contract NAS5-26555. PB was partially funded by a Packard Fellowship, Grant #2019-69646. PB was also partially supported by a Giacconi Fellowship from the Space Telescope Science Institute. Finally, PB was also partially supported through program number HST-HF2-51353.001-A, provided by NASA through a Hubble Fellowship grant from the Space Telescope Science Institute, under NASA contract NAS5-26555.

Data compilations from many studies used in this paper were made much more accurate and efficient by the online WEBPLOTDIGITIZER code.<fnref n="3" /> This research has made extensive use of the arXiv and NASA's Astrophysics Data System.

This research used the Ocelote supercomputer of the University of Arizona. The allocation of computer time from the UA Research Computing High Performance Computing (HPC) at the University of Arizona is gratefully acknowledged. The Bolshoi-Planck simulation was performed by Anatoly Klypin within the Bolshoi project of the University of California High-Performance AstroComputing Center (UC-HiPACC; PI Joel Primack).

# DATA AVAILABILITY

The parallel implementation of TRINITY, the compiled data sets (Section 3.2), the posterior distribution of model parameters, and the codes to compute the predictions for all the plots in this work are available at https://github.com/HaowenZhang/TRINITY.

<footnote n="3">³https://apps.automeris.io/wpd/</footnote>
REFERENCES

Aird J. et al., 2010, MNRAS, 401, 2531

Aird J., Coil A. L., Georgakakis A., 2018, MNRAS, 474, 1225

Alexander D. M., Hickox R. C., 2012, New Astron. Rev., 56, 93

Anglés-Alcázar D., Faucher-Giguère C.-A., Quataert E., Hopkins P. F., Feldmann R., Torrey P., Wetzel A., Kereš D., 2017, MNRAS, 472, L109

Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143

Bouwens R. J., Stefanon M., Oesch P. A., Illingworth G. D., Nanayakkara T., Roberts-Borsani G., Labbé I., Smit R., 2019, ApJ, 880, 25

Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

Bower R. G., Schaye J., Frenk C. S., Theuns T., Schaller M., Crain R. A., McAlpine S., 2017, MNRAS, 465, 32

Brandt W. N., Alexander D. M., 2015, A&AR, 23, 1

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Bryan G. L., Norman M. L., 1998, ApJ, 495, 80

Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storch-Bergmann T., 2000, ApJ, 533, 682

Caplar N., Lilly S. J., Trakhtenbrot B., 2015, ApJ, 811, 148

Caplar N., Lilly S. J., Trakhtenbrot B., 2018, ApJ, 867, 148

Cavagnolo K. W., McNamara B. R., Nulsen P. E. J., Carilli C. L., Jones C., Bîrzan L., 2010, ApJ, 720, 1066

Chabrier G., 2003, PASP, 115, 763

Choi E., Ostriker J. P., Naab T., Somerville R. S., Hirschmann M., Núñez A., Hu C.-Y., Oser L., 2017, ApJ, 844, 31

Cirasuolo M., Magliocchetti M., Celotti A., 2005, MNRAS, 357, 1267

Conroy C., White M., 2013, ApJ, 762, 70

Costa T., Sijacki D., Haehnelt M. G., 2014, MNRAS, 444, 2355

Croton D. J. et al., 2006, MNRAS, 365, 11

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827

Delvecchio I. et al., 2014, MNRAS, 439, 2736

Dubois Y., Volonteri M., Silk J., 2014, MNRAS, 440, 1590

Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., Teyssier R., 2015, MNRAS, 452, 1502

Fanidakis N. et al., 2012, MNRAS, 419, 2797

Ferraese L., Merritt D., 2000, ApJ, 539, L9

Gebhardt K. et al., 2000, ApJ, 539, L13

Gültekin K. et al., 2009, ApJ, 698, 198

Haario H., Saksman E., Tamminen J., 2001, Bernoulli, 7, 223

Habouzit M., Volonteri M., Dubois Y., 2017, MNRAS, 468, 3935

Habouzit M. et al., 2022, MNRAS, 509, 3015

Harikane Y., Nakajima K., Ouchi M., Umeda H., Isobe Y., Ono Y., Xu Y., Zhang Y., 2024, ApJ, 960, 56

Häring N., Rix H.-W., 2004, ApJ, 604, L89

Heckman T. M., Best P. N., 2014, ARA&A, 52, 589

Hirschmann M., Somerville R. S., Naab T., Burkert A., 2012, MNRAS, 426, 237

Hlavacek-Larrondo J. et al., 2015, ApJ, 805, 35

Ho L. C., 2002, ApJ, 564, 120

Ho L. C., 2008, ARA&A, 46, 475

Hopkins P. F., Bundy K., Hernquist L., Ellis R. S., 2007, ApJ, 659, 976

Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2018, ApJ, 854, 73

Izumi T. et al., 2019, PASJ, 71, 111

Jarosik N. et al., 2011, ApJS, 192, 14

Kelly B. C., Shen Y., 2013, ApJ, 764, 45

Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102

Kormendy J., Ho L. C., 2013, ARA&A, 51, 511

Kormendy J., Richstone D., 1995, ARA&A, 33, 581

La Franca F., Melini G., Fiore F., 2010, ApJ, 718, 368

Lang P. et al., 2014, ApJ, 788, 11

Lauer T. R., Tremaine S., Richstone D., Faber S. M., 2007, ApJ, 670, 249

Magorrian J. et al., 1998, AJ, 115, 2285

Maiolino R. et al., 2023, Nature, preprint (arXiv:2305.12492)

Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, MNRAS, 351, 169

McConnell N. J., Ma C.-P., 2013, ApJ, 764, 184

McDonald M., McNamara B. R., Calzadilla M. S., Chen C.-T., Gaspari M., Hickox R. C., Kara E., Korchagin I., 2021, ApJ, 908, 85

Mendel J. T., Simard L., Palmer M., Ellison S. L., Patton D. R., 2014, ApJS, 210, 3

Merloni A., Heinz S., 2007, MNRAS, 381, 589

Merloni A., Heinz S., 2008, MNRAS, 388, 1011

Merloni A., Rudnick G., Di Matteo T., 2004, MNRAS, 354, L37

Mineshige S., Kawaguchi T., Takeuchi M., Hayashida K., 2000, PASJ, 52, 499

Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121

Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822

Nagar N. M., Falcke H., Wilson A. S., 2005, A&A, 435, 521

Narayan R., Yi I., 1994, ApJ, 428, L13

Oesch P. A., Bouwens R. J., Illingworth G. D., Labbé I., Stefanon M., 2018, ApJ, 855, 105

Peca A. et al., 2023, ApJ, 943, 162

Planck Collaboration XIII, 2016, A&A, 594, A13

Reines A. E., Volonteri M., 2015, ApJ, 813, 82

Rodríguez-Puebla A., Behroozi P., Primack J., Klypin A., Lee C., Hellinger D., 2016, MNRAS, 462, 893

Savorgnan G. A. D., Graham A. W., Marconi A. R., Sani E., 2016, ApJ, 817, 21

Schulze A., Wisotzki L., 2010, A&A, 516, A87

Schulze A. et al., 2015, MNRAS, 447, 2085

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shankar F., Cavaliere A., Cirasuolo M., Maraschi L., 2008, ApJ, 676, 131

Shankar F., Weinberg D. H., Miralda-Escudé J., 2009, ApJ, 690, 20

Shankar F., Weinberg D. H., Miralda-Escudé J., 2013, MNRAS, 428, 421

Shankar F. et al., 2020, MNRAS, 493, 1500

Shen Y., 2013, Bull. Astron. Soc. India, 41, 61

Shen X., Hopkins P. F., Faucher-Giguère C.-A., Alexander D. M., Richards G. T., Ross N. P., Hickox R. C., 2020, MNRAS, 495, 3252

Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877

Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575

Silk J., Rees M. J., 1998, A&A, 331, L1

Silverman J. D. et al., 2008, ApJ, 679, 118

Smolčić V. et al., 2017, A&A, 602, A6

Soltan A., 1982, MNRAS, 200, 115

Somerville R. S., Hopkins P. F., Cox T. J., Robertson B. E., Hernquist L., 2008, MNRAS, 391, 481

Tillman M. T., Wellons S., Faucher-Giguère C.-A., Kelley L. Z., Anglés-Alcázar D., 2022, MNRAS, 511, 5756

Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709

Tremaine S. et al., 2002, ApJ, 574, 740

Tremmel M., 2017, PhD thesis, University of Washington, Seattle

Ueda Y., Akiyama M., Hasinger G., Miyaji T., Watson M. G., 2014, ApJ, 786, 104

Vogelsberger M. et al., 2014, MNRAS, 444, 1518

Volonteri M., Dubois Y., Pichon C., Devriendt J., 2016, MNRAS, 460, 2979

Weinberger R. et al., 2017, MNRAS, 465, 3291

Yang G. et al., 2018, MNRAS, 475, 1887

Yang J. et al., 2021, ApJ, 923, 262

Yu Q., Tremaine S., 2002, MNRAS, 335, 965

Zhang H., Behroozi P., Volonteri M., Silk J., Fan X., Hopkins P. F., Yang J., Aird J., 2023a, MNRAS, 523, L69

Zhang H. et al., 2023b, MNRAS, preprint (arXiv:2309.07210)

Zhang H., Behroozi P., Volonteri M., Silk J., Fan X., Hopkins P. F., Yang J., Aird J., 2023c, MNRAS, 518, 2123
# APPENDIX: RADIATIVE AND TOTAL EFFICIENCIES

## A1 The fiducial scaling relation between $\epsilon_{rad}$ and $\epsilon_{tot}$

In TRINITY, we adopt a non-linear scaling relation between radiative and kinetic components of AGN energy efficiency, $\epsilon_{rad}$ and $\epsilon_{kin}$, and the total energy efficiency, $\epsilon_{tot}$. This is motivated by the fact that active SMBHs convert gravitational energy not only into radiation, but also into kinetic energy in the form of jets and/or outflows. The scaling relation is as follows:

$$ \epsilon_{\text{rad}} = \epsilon_{\text{tot}} \times \begin{cases} \frac{\eta}{\eta_{\text{crit}}}, & \eta \le \eta_{\text{crit}} \\ \frac{1}{2/\eta[1 + \ln(\eta/2)]}, & \eta_{\text{crit}} < \eta \le 2, \\ \eta > 2 & \eta > 2 \end{cases} \quad (\text{A1}) $$

and:

$$ \epsilon_{\text{kin}} = \begin{cases} \epsilon_{\text{tot}}(1 - \eta/0.03), & \eta \le \eta_{\text{crit}} \\ 0, & \eta > \eta_{\text{crit}} \end{cases}, \quad (\text{A2}) $$

where $\eta_{crit} = 0.03$. This is also equivalent to defining radiative and kinetic components of Eddington ratio, $\eta_{rad}$ and $\eta_{kin}$, so that $\eta = \eta_{rad} + \eta_{kin}$:

$$ \eta_{\text{rad}} = \begin{cases} \frac{\eta^2}{\eta_{\text{crit}}}, & \eta \le \eta_{\text{crit}} \\ \eta, & \eta_{\text{crit}} < \eta \le 2, \\ 2[1 + \ln(\eta/2)], & \eta > 2 \end{cases}, \quad (\text{A3}) $$

and

$$ \eta_{\text{kin}} = \begin{cases} \eta - \frac{\eta^2}{\eta_{\text{crit}}}, & \eta \le \eta_{\text{crit}} \\ 0, & \eta > \eta_{\text{crit}} \end{cases}. \quad (\text{A4}) $$

As is shown in Fig. A1, the scaling between $\eta_{rad}$ and $\eta$ is similar to the one used by Merloni & Heinz (2008) at $\eta \le 2$, which was based on the observed scaling found by Merloni & Heinz (2007). But Merloni & Heinz (2008) adopted multiple scaling relations between $\eta$, $\eta_{rad}$, and $\eta_{kin}$, each describing a different AGN accretion and feedback mode. At $\eta > \eta_{crit}$ the “high kinetic” (HK, solid lines) mode features similar kinetic energy and radiation outputs, and radiation dominates over kinetic energy in the “high radiative” mode (HR, dashed lines). Merloni & Heinz (2008) further assumed that the ratio of HK and HR quasars is 1:9. The $\eta_{rad}-\eta$ relations were derived with observed scaling relations between AGN radio and X-ray luminosities, and SMBH masses, which came with substantial scatter. In light of this, we choose to stick to the simpler scaling in equation (A3). For $\eta \ge 2$, we adopt a logarithmic function to account for the part of the outgoing radiation is trapped by the thickened accretion disc (Mineshige et al. 2000). At a given $\eta < \eta_{crit}$ (“low kinetic” mode, or LK in Merloni & Heinz 2008), equation (A4) produces more kinetic energy than Merloni & Heinz (2008). We also ignore the kinetic energy output from active SMBHs with $\eta > \eta_{crit}$, due to the lack of observational constraints.

## A2 Variants of the $\epsilon_{rad}-\epsilon_{tot}$ relation

To account for the systematic uncertainties induced by our assumed $\epsilon_{rad}-\epsilon_{tot}$ relation, we also reran TRINITY with two different variants of the $\epsilon_{rad}-\epsilon_{tot}$ relation: (1) the one where $\epsilon_{rad} \equiv \epsilon_{tot}$, i.e. no kinetic energy is produced from SMBH accretion (“no kinetic” hereafter), and (2) the one where $\eta_{crit} = 0.1$, i.e. more kinetic energy is generated in SMBH accretion (“more kinetic” hereafter). These two variants represent the two extremes of kinetic energy output in SMBH accretion. All the qualitative conclusions remain unchanged with both variant models. To save space, we uploaded Figs 2–12 from these two variant models https://github.com/HaowenZhang/TRINITY/tree/main/plots/Paper3/Model_Variants.

<figure><img src="image_17.png" /><figcaption>Figure A1. The comparison between the $\eta_{rad}-\eta_{kin}-\eta_{tot}$ scaling relations used in TRINITY and Merloni & Heinz (2008). The orange solid line is the 1:1 equality line. Three AGN accretion modes are defined by Merloni & Heinz (2008): the “low kinetic” (LK, solid lines) mode with $\eta_{tot} \le 0.03$, the “high kinetic” (HK, solid lines) mode, and the “high radiative” (HR, dashed lines) mode. An AGN with $\eta_{tot} > 0.03$ is in either the HK or HR mode. $\eta_{rad}$ and $\eta_{kin}$ are shown in blue and red curves, respectively. The scaling relations used in TRINITY are shown in dotted lines.</figcaption></figure>

This paper has been typeset from a TeX/LATEX file prepared by the author.