# TRINITY IV: predictions for supermassive black holes at z 6

Haowen Zhang, Peter Behroozi, Marta Volonteri, Joseph Silk, Xiaohui Fan, James Aird, Jinyi Yang, Feige Wang, Wei Leong Tee, Philip F. Hopkins

► To cite this version:

Haowen Zhang, Peter Behroozi, Marta Volonteri, Joseph Silk, Xiaohui Fan, et al. TRINITY IV: predictions for supermassive black holes at z 6. Monthly Notices of the Royal Astronomical Society, 2024, 531, pp.4974-4989. 10.1093/mnras/stae1447. insu-04851547

HAL Id: insu-04851547

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

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
# TRINITY IV: predictions for supermassive black holes at $z \gtrsim 6$

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="1" />,<fnref n="3" />,<fnref n="4" />,<fnref n="5" /> Xiaohui Fan,<fnref n="1" />

James Aird<fnref n="1" />,<fnref n="6" />,<fnref n="7" /> Jinyi Yang (杨锦怡),<fnref n="1" /> Feige Wang (王飞格),<fnref n="1" /> Wei Leong Tee<fnref n="1" /> and Philip F. Hopkins<fnref n="1" />,<fnref n="8" />

Accepted 2024 May 31. Received 2024 May 31; in original form 2023 September 5

## ABSTRACT

We present predictions for the high-redshift halo–galaxy–supermassive black hole (SMBH) connection from the TRINITY model. Matching a comprehensive compilation of galaxy ($0 \le z \le 13$) and SMBH data sets ($0 \le z \le 6.5$), TRINITY finds: (1) The number of SMBHs with $M_{\bullet} > 10^9$ M$_{\odot}$ in the observable Universe increases by five orders of magnitude from $z \sim 10$ to $z \sim 2$, and by another factor of $\sim 3$ from $z \sim 2$ to $z = 0$; (2) The $M_{\bullet} > 10^9$ and $10^{10}$ M$_{\odot}$ SMBHs at $z \sim 6$ live in haloes with $\sim (2-3)$ and $(3-5) \times 10^{12}$ M$_{\odot}$; (3) the newly discovered JWST AGN candidates at $7 \le z \le 11$ are overmassive compared to the intrinsic SMBH mass–galaxy mass relation from TRINITY, but they are still broadly consistent with TRINITY predictions for flux limited AGN samples with Lauer bias. This bias favours the detection for overmassive SMBHs due to higher luminosities at a fixed Eddington ratio. However UHZ1’s $M_{\bullet}/M_{\odot}$ ratio is still some 1 dex higher than TRINITY AGNs, indicating a discrepancy; (4) TRINITY underpredicts the number densities of GN-z11 and CEERS_1019 analogues. But given the strong constraints from existing data in TRINITY, the extra constraint from GN-z11 and CEERS_1019 does not significantly change TRINITY model results. ($z = 6–10$ quasar luminosity functions will reduce uncertainties in the TRINITY prediction of the $z = 6–10$ SMBH mass–galaxy mass relation by up to $\sim 0.5$ dex. These luminosity functions will be available with future telescopes, such as Roman and Euclid.

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

## 1 INTRODUCTION

Supermassive black holes (SMBHs) are believed to exist in the centres of most, if not all 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 SMBHs are actively accreting, they are called active galactic nuclei (AGN). AGNs are also called quasars when they outshine their host galaxies. With their accretion energy output, SMBHs are likely able to modulate both their own mass accretion and the star formation of their host galaxies. (Silk & Rees 1998; Bower et al. 2006; Somerville et al. 2008; Sijacki et al. 2015). On the other hand, galaxies can also regulate SMBH growth via the fueling of gas and galaxy mergers. Consequently, both SMBHs and host galaxies may influence each others’ growth, which is known as ‘coevolution’. Hence, to understand galaxy and SMBH assembly histories, it is critical to constrain the interaction between SMBHs and their host galaxies (see e.g. Hopkins et al.

2007; Ho 2008; Alexander & Hickox 2012; Kormendy & Ho 2013; Heckman & Best 2014; Brandt & Alexander 2015).

Many observations support the coevolution scenario. For example, there are relatively tight scaling relations ($\sim 0.3$ dex scatter) between SMBH masses and host galaxy dynamical properties (e.g. velocity dispersion or bulge mass) at $z \sim 0$ (see Häring & Rix 2004; Gültekin et al. 2009; Kormendy & Ho 2013; McConnell & Ma 2013; Savorgnan et al. 2016). In addition, the cosmic SMBH accretion rate (CBHAR) density is roughly proportional to cosmic star formation rate (CSFR) density over $0 < z < 4$, with a ratio of $\sim 10^{-4}-10^{-3}$ (Merloni 2004; Silverman et al. 2008; Shankar, Weinberg & Miralda-Escudé 2009; Aird et al. 2010; Delvecchio et al. 2014; Yang et al. 2018). However, it has been difficult to verify other predictions of the coevolution model (e.g. a tight SMBH–galaxy mass relation) at higher redshifts.

At $z = 0$, high spatial resolution spectroscopy and dynamics modelling have been used to measure SMBH–galaxy scaling relations (e.g. Magorrian et al. 1998; Ferrarese & Ford 2005; McConnell & Ma 2013). But it is difficult to measure individual SMBH masses in the same way at $z > 0$ due to lower spatial resolutions. Hence, we have to rely on indirect methods such as reverberation mapping (Blandford & McKee 1982; Peterson 1993) or ‘virial estimates’ with

<footnote n="1">¹Department of Astronomy, University of Arizona, 933 N Cherry Ave., 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 Université), 98 bis Bd. Arago, F-75014 Paris, 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>
single-epoch AGN spectra, which infer SMBH masses from AGN emission line width and luminosity (Vestergaard & Peterson 2006). All such indirect methods require that SMBHs are actively accreting, and thus: (1) are biased against dormant SMBHs, and (2) make it difficult to measure host galaxy masses at the same time. In addition, overmassive SMBHs are brighter and more likely to be included in the sample, when different SMBHs have similar Eddington ratios. This bias affects the measurement of the galaxy mass–SMBH mass ($M_{\bullet}-M_{*}$) relation with bright quasars (Lauer et al. 2007). Zhang et al. (2023a) quantified this bias as a function of AGN luminosity, and showed that AGNs with bolometric luminosities $L_{bol} \lesssim 10^{45}$ erg s$^{-1}$ are needed to measure the unbiased $M_{\bullet}-M_{*}$ relation at $z = 6$.

Prior to JWST, the highest-redshift AGN was found at $z \sim 7.6$ (Wang et al. 2021). AGN luminosity functions at up to $z \sim 7$ have also been measured (Kulkarni, Worseck & Hennawi 2019; Wang et al. 2019; Shen et al. 2020). Most $z \sim 7$ AGNs were found to be bright quasars with high-SMBH masses, accreting at around the Eddington rate (Shen et al. 2019). However, our observing capabilities made it difficult to ascertain if such mass and Eddington ratio distributions were representative of all high-redshift SMBHs.

JWST has changed the landscape of high-redshift quasar observations. As of the writing of this paper, JWST has found at least two AGN candidates above $z \sim 8$, i.e. CEERS_1019 at $z = 8.7$ (Larson et al. 2023) and GN-z11 at $z = 10.6$ (Maiolino et al. 2023a). CEERS_1019 has an estimated SMBH mass $M_{\bullet} \sim 10^7 M_{\odot}$, and an Eddington ratio $\eta \sim 1.3$, whereas GN-z11 is estimated to be a $M_{\bullet} \sim 10^6 M_{\odot}$ SMBH accreting at an Eddington ratio of $\eta \sim 5$. JWST's sensitivity finally allows measuring SMBH masses in AGNs with $L_{bol} \lesssim 10^{45}$ erg s$^{-1}$ at $z \gtrsim 6$ (e.g. CEERS_1019 and GN-z11), as well as measuring host galaxy properties for high-z X-ray AGN candidates, like UHZ-1 at $z = 10.1$ (Bogdan et al. 2023; Goulding et al. 2023). Therefore, JWST will continue to provide strong tests/constraints on theoretical models of the halo–galaxy–SMBH connection.

Different theoretical studies tried to ascertain the possible origins and evolution histories of these highest-redshift AGNs. Using a semi-analytical model, the Cosmic Archaeology Tool (CAT; Trinca et al. 2022), Schneider et al. (2023) found a potential evolutionary connection between objects like GN-z11 and CEERS_1019. Scholtz et al. (2023) estimated GN-z11's host halo mass to be $\sim 3 \times 10^{10} M_{\odot}$, based on the stellar mass–halo mass relation from the empirical UNIVERSE MACHINE model (Behroozi et al. 2019). This halo mass scale makes GN-z11 a candidate protocluster core, which could reach $\gtrsim 10^{15} M_{\odot}$ at $z = 0$.

In this paper, we present the SMBH properties at $z \gtrsim 6$ from the empirical TRINITY model. We further put CEERS_1019 and GN-z11 in the context of TRINITY, and examine their (in)consistency with our predictions. TRINITY makes statistical connections between dark matter haloes, galaxies, and SMBHs from $z = 0-10$. Specifically, TRINITY infers SMBH growth histories as a function of halo/galaxy mass in a similar manner as applying the Soltan Argument (Soltan 1982) to each halo/galaxy population (see also Shankar et al. 2020). Compared to previous empirical models, TRINITY includes joint constraints from a large and diverse compilation of data. Specifically, the following data are used to constrain the halo–galaxy connection: 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). Additionally, the galaxy–SMBH connection is constrained by quasar luminosity functions (QLFs), quasar probability distribution functions (QPDFs), active black hole mass functions (ABHMFs), and the $z = 0$ bulge mass–black hole mass relations. To accommodate the enormous joint constraining power of this data set, TRINITY has a more flexible

parametrization than past approaches (e.g. Marconi et al. 2004; Shankar et al. 2009, 2013; Conroy & White 2013). Lastly, TRINITY features more realistic modeling of AGN observables by including, e.g. black hole mergers and kinetic AGN luminosities in the model.

Following Behroozi, Wechsler & Conroy (2013), TRINITY modelling starts with halo population statistics from an N-body simulation. The model then generates mock universes, each with different galaxy and SMBH mass–growth rate distributions as functions of halo mass and redshift. Millions of mock universes are compared with real observations using a Markov Chain Monte Carlo (MCMC) algorithm to converge on the range of halo–galaxy–SMBH connections that matches all the input data.

The paper is organized as follows. In Section 2, we briefly introduce the TRINITY framework. Section 3 covers the N-body simulations and galaxy/SMBH connection. For full details about the TRINITY methodology, simulations, and data sets, we refer readers to Zhang et al. (2023c). In Section 4, we present high-redshift predictions from TRINITY. Section 5 compares CEERS1019 and GN-z11 with TRINITY's predictions, and discusses how future observations would help constrain the demographics and growth histories of high redshift SMBHs. Finally, our conclusions are given 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 the key aspects of the TRINITY high-redshift halo–galaxy–SMBH connection that induce the predictions in this work, and how they are constrained by our data compilation. We then give a brief overview on TRINITY's implementation in Section 2.2. For full details, we refer readers to Section 2 of Zhang et al. (2023c).

### 2.1 Inferring the high-redshift halo–galaxy–SMBH connection with TRINITY

TRINITY infers the range of empirical halo–galaxy–SMBH connections that match all of our galaxy+AGN data. Many predictions in this work come from the following two aspects of this connection:

First, the intrinsic $M_{\bullet}-M_{*}$ relation evolves only mildly with redshift (see Section 4.1). The largest change occurs for low-mass galaxies (i.e. $M_{*} \lesssim 5 \times 10^9 M_{\odot}$) which have lower SMBH masses at higher redshifts; almost no change in SMBH masses with redshift occurs for galaxies with masses $M_{*} \gtrsim 5 \times 10^{10} M_{\odot}$. This redshift evolution is mainly constrained by the combination of QLFs, QPDFs, and the local $M_{\bullet}-M_{\text{bulge}}$ relation. Specifically, different redshift evolutions of the $M_{\bullet}-M_{*}$ relation produce different average SMBH growth rates as functions of host galaxy mass and redshift. Therefore, there is only one unique combination of a redshift-dependent $M_{\bullet}-M_{*}$ relation and an AGN energy efficiency to simultaneously match: (1) the AGN luminosity distribution as a function of galaxy mass and redshift, i.e. QPDFs; (2) QLFs; and (3) the local $M_{\bullet}-M_{\text{bulge}}$ relation. All other redshift-dependent $M_{\bullet}-M_{*}$ relations will cause tension with at least some data constraints. For example, a steeper $M_{\bullet}-M_{*}$ relation will violate the constraint from the local $M_{\bullet}-M_{\text{bulge}}$ relation, induce
faster growth of SMBHs at fixed galaxy mass, and overproduce AGN activity compared to observed QPDFs.

Secondly, the AGN Eddington ratio distribution becomes narrower towards higher redshifts, with increasing average Eddington ratios. This is required by the ABHMFs from Kelly & Shen (2013). With the redshift-dependent $M_{\bullet}-M_{*}$ relation and AGN duty cycles constrained by QPDFs, QLFs, and the local $M_{\bullet}-M_{\text{bulge}}$ relation, lower Eddington ratios and/or broader Eddington ratio distributions will underpredict the number of SMBHs above certain luminosity limits, i.e. ABHMFs.

Of note, the AGN data for constraining TRINITY spans a redshift range of $z = 0-6.5$. When making predictions for SMBHs beyond $z \sim 7$, we are effectively examining the observational consequences of a mildly evolving SMBH–galaxy relationship at those redshifts. Hence, any inconsistency between TRINITY predictions with $z > 7$ observations would indicate different evolution in the SMBH–galaxy mass relationship than that seen at $z < 7$. As we will show in Sections 4.1 and 5.2, the SMBH properties of CEERS_1019 and GN-z11 are partly consistent with these extrapolations. The partial inconsistency between TRINITY and these two AGNs does not impose strong constraints on the current TRINITY model, i.e. the inconsistency is not yet statistically significant. As a result, further observations of high-redshift AGNs are needed to either confirm our prediction, or otherwise provide stronger constraints on the TRINITY model.

## 2.2 Implementation overview

TRINITY makes self-consistent and statistical halo–galaxy–SMBH connections from $z = 0-10$. Instead of following individual haloes across cosmic time like the UNIVERSEMACHINE (Behroozi et al. 2019), TRINITY is built upon halo statistics from N-body simulations. The halo–galaxy connection is made via the SFR–$M_{\text{peak}}$ relation, where $M_{\text{peak}}$ is the peak halo mass of haloes throughout their assembly histories. With a given redshift-dependent SFR–$M_{\text{peak}}$ relation, TRINITY calculates the average star formation histories for different halo populations, and integrates them over time to obtain average galaxy stellar masses ($M_{*}$) as functions of redshift and halo mass. Average galaxy total stellar masses are then converted into average bulge masses ($M_{\text{bulge}}$) with a redshift-dependent but fixed $M_{\text{bulge}}-M_{*}$ based on SDSS (Mendel et al. 2014) and CANDELS (Lang et al. 2014) observations.

To connect SMBHs with host galaxies, We parametrize the $M_{\bullet}-M_{\text{bulge}}$ relation as a redshift-dependent power-law for all galaxies, i.e. SMBH hosts and non-hosts. The lognormal scatter around the $M_{\bullet}-M_{\text{bulge}}$ relation, $\sigma_{\text{BH}}$, is a free parameter to be constrained in MCMC. We assume that $\sigma_{\text{BH}}$ is redshift-independent due to the lack of data to constrain its redshift dependence. As shown in Zhang et al. (2023c), our results do not change significantly if we instead parametrize the $M_{\bullet}-M_{*}$ relation as a power law. SMBH occupation fractions are parametrized as a sigmoid function of halo mass, which evolves with redshift. Together, these choices fully parametrize the SMBH mass distribution as a function of galaxy mass and redshift ($P(M_{\bullet}|M_{*}, z)$). Average black hole accretion rates (BHARs) are calculated by comparing the average $M_{\bullet}$ of the same halo population between two consecutive snapshots.

To forward model SMBH observables, we convert average BHARs into AGN Eddington ratio distributions using parametrized AGN efficiency, duty cycles, and Eddington ratio distribution shapes (i.e. $P(L_{\text{bol}}|M_{\bullet}, M_{*}, z)$). Before comparing predicted observables with real data, we also correct them for systematic effects like errors in stellar masses and SFRs from SED fitting, and Compton-thick AGN obscuration. To get the posterior distribution of model parameters, we compare predicted observables with the compiled data for $\sim 2$

million points in parameter space with a custom made Metropolis Markov Chain Monte Carlo (MCMC) algorithm (based on Haario, Saksman & Tamminen 2001). From this posterior distribution, we obtain both the best-fitting model parameter set given the observational constraints, as well as the uncertainties in the halo–galaxy–SMBH connection and AGN properties, e.g. efficiency and duty cycles. These uncertainties quantify how the predicted halo–galaxy–SMBH connection and SMBH properties are can vary without violating our input data constraints. The whole parametrization of TRINITY contains 54 free parameters.

## 2.3 Calculating SMBH statistics

In this work, we present various $z \gtrsim 7$ SMBH statistics from the best-fitting TRINITY model. Here, we briefly introduce how we calculate them.

To calculate the median $M_{\bullet}-M_{*}$ relation for bright quasars above a certain bolometric luminosity threshold ($L_{\text{bol, lim}}$), we first calculate the SMBH mass function at each given stellar mass and redshift:

$$ \phi_{\bullet}(M_{\bullet}|M_{*}, L_{\text{bol}} > L_{\text{bol,lim}}, z) = \int_{0}^{\infty} \phi_{\text{h}}(M_{\text{peak}}, z) \int_{L_{\text{bol,lim}}}^{\infty} P(M_{\bullet}|M_{*}, z)P(M_{*}|M_{\text{peak}}, z) P(L_{\text{bol}}|M_{\text{peak}}, z) dL_{\text{bol}}dM_{\text{peak}}, \quad (1) $$

where $P(L_{\text{bol}}|M_{*}, z)$ is obtained by integrating the number density of black holes emitting at the corresponding Eddington ratio:

$$ P(L_{\text{bol}}|M_{\text{peak}}, z) = \int_{0}^{\infty} P(\eta_{\text{rad}}(L_{\text{bol}}, M_{\bullet}) | M_{\bullet}, M_{\text{peak}}, z) \times P(M_{\bullet}|M_{\text{peak}}, z) dM_{\bullet}. \quad (2) $$

We then calculate the median $M_{\bullet}$ and the 16th–84th percentile range from $\phi_{\bullet}(M_{\bullet}|M_{*}, L_{\text{bol}} > L_{\text{bol, lim}}, z)$. The $M_{\bullet}-M_{*}$ relation for bright quasars is then obtained by repeating this process for a series of $M_{*}$ at a fixed $z$.

$\phi_{\bullet}(M_{\bullet}|M_{*}, L_{\text{bol}} > L_{\text{bol, lim}}, z)$ can also be used for calculating the expected number of AGNs above $L_{\text{bol, lim}}$ and observed SMBH mass $M_{\bullet, \text{obs}}$ in galaxies with observed galaxy mass $M_{*\text{, obs}}$, galaxy mass uncertainty $\sigma_{M_{*}}$, and SMBH mass uncertainty $\sigma_{M_{*}}$. We first calculate the cumulative AGN number density above $L_{\text{bol, lim}}$ and $M_{\bullet, \text{obs}}$, $\phi_{\bullet}$, as a function of $M_{*\text{, obs}}$, $\sigma_{M_{*}}$, and redshift:

$$ \phi_{>} = \phi_{>} (M_{*\text{, obs}}, \sigma_{M*}, z) = \int_{0}^{\infty} P(M_{\bullet} > M_{\bullet, \text{obs}} | M_{\bullet, \text{obs}}, \sigma_{M*}) dM_{\bullet} \int_{0}^{\infty} \phi_{\bullet}(M_{\bullet}|M_{*}, L_{\text{bol}} > L_{\text{bol,lim}}, z) \times P(M_{*}|M_{*\text{, obs}}, \sigma_{M*}) dM_{*}, \quad (3) $$

where $P(M_{*} > M_{*\text{, obs}}|M_{*\text{, obs}}, \sigma_{M*})$ is the cumulative probability that the intrinsic SMBH mass is larger than the observed value. Assuming a lognormal random scatter of $\sigma_{M*} = 0.5$ dex (see e.g. McLure & Dunlop 2004 and Vestergaard & Osmer 2009), $P(M_{*} > M_{*\text{, obs}}|M_{*\text{, obs}}, \sigma_{M*})$ is a complementary error function centered at $M_{*\text{, obs}}$. $P(M_{*}|M_{*\text{, obs}}, \sigma_{M*})$ is a lognormal distribution centered on $M_{*\text{, obs}}$ with a spread of $\sigma_{M*}$. Between $z_1 < z < z_2$, the expected number of AGNs above $L_{\text{bol, lim}}$ and $M_{*\text{, obs}}$ in galaxies with $M_{*\text{, obs}}$ and $\sigma_{M*}$ is:

$$ N_{\text{exp}}(M_{*\text{, obs}}, \sigma_{M*}, z_1 < z < z_2) = \int_{z_1}^{z_2} \phi_{>} (M_{*\text{, obs}}, \sigma_{M*}, z) \frac{dV}{dz} dz, \quad (4) $$
where $dV/dz$ is the comoving volume between $z_1$ and $z_2$ over the whole sky. For any given survey, $dV/dz$ should be further corrected with the actual survey area and any lensing magnification, if applicable. Assuming Poisson statistics, the probability of detecting *at least* one such AGN is:

$$P(N_{\text{obs}} \ge 1|N_{\text{exp}}) = 1 - P(N_{\text{obs}} = 0|N_{\text{exp}}) = 1 - e^{-N_{\text{exp}}}. \quad (5)$$

SMBH mass functions are convolutions of halo mass functions and the conditional probability distribution of $M_\bullet$ given $M_{\text{peak}}$:

$$\phi(M_{\bullet}, z) = \int_{0}^{\infty} \phi_{h}(M_{\text{peak}}, z) P(M_{\bullet}|M_{\text{peak}}, z) dM_{\text{peak}}, \quad (6)$$

where $P(M_\bullet|M_{\text{peak}}, z)$ is the convolution of: (1) the probability distribution of $M_\bullet$ given $M_*$, $P(M_\bullet|M_*, z)$; and (2) the probability distribution of $M_\bullet$ given $M_{\text{peak}}$, $P(M_\bullet|M_{\text{peak}}, z)$:

$$P(M_{\bullet}|M_{\text{peak}}, z) = \int_{0}^{\infty} P(M_{\bullet}|M_{*}, z)P(M_{*}|M_{\text{peak}}, z)dM_{*}. \quad (7)$$

The cosmic number density of SMBHs above a mass threshold $M_{\bullet, \text{lim}}$, $n_\bullet(M_\bullet > M_{\bullet, \text{lim}}, z)$, is a convolution of $P(M_\bullet|M_{\text{peak}}, z)$ and the halo mass function $\phi_h(M_{\text{peak}})$:

$$n_{\bullet}(M_{\bullet} > M_{\bullet, \text{lim}}, z) = \int_{M_{\bullet, \text{lim}}}^{\infty} P(M_{\bullet}|M_{\text{peak}}, z) \phi_h(M_{\text{peak}}, z) dM_{\text{peak}}. \quad (8)$$

The differential number of SMBHs above $M_{\bullet, \text{lim}}$ per unit redshift, $dN_\bullet(M_\bullet > M_{\bullet, \text{lim}}, z)/dz$, is then:

$$\frac{dN_{\bullet}}{dz}(M_{\bullet} > M_{\bullet, \text{lim}}, z) = n_{\bullet}(M_{\bullet} > M_{\bullet, \text{lim}}, z) \times \frac{dV}{dz}, \quad (9)$$

where $dV/dz$ is the differential comoving volume per unit redshift, which is calculated based on the adopted cosmology (see Section 1). The cumulative number of SMBHs above $M_{\bullet, \text{lim}}$ at above redshift $z$, $N_\bullet(M_\bullet > M_{\bullet, \text{lim}}, >z)$ is thus:

$$N_{\bullet}(M_{\bullet} > M_{\bullet, \text{lim}}, > z) = \int_{z}^{\infty} \frac{dN_{\bullet}}{dz}(M_{\bullet} > M_{\bullet, \text{lim}}, z')dz'. \quad (10)$$

# 3 SIMULATIONS AND DATA CONSTRAINTS

## 3.1 Dark matter halo statistics

As noted in Section 2.2, TRINITY requires only halo population statistics from dark matter simulations, as opposed to individual halo merger trees. We use the peak historical mass ($M_{\text{peak}}$) halo mass functions for central and satellite haloes from Behroozi et al. (2013), for the cosmology specified in the introduction. These mass functions are based on central halo mass functions from Tinker et al. (2008), with adjustments to include satellite halo number densities as well as to use $M_{\text{peak}}$ instead of the present day halo mass. These adjustments were based on the Bolshoi & Consuelo simulations (Klypin, Trujillo-Gomez & Primack 2011). We refer readers to appendix G of Behroozi et al. (2013) for full details. Note that, the WMAP7 cosmology (Jarosik et al. 2011) was used in the Bolshoi simulation, whereas TRINITY assumes the Planck cosmology. We thus made further corrections to halo mass functions for this cosmology difference. We calculate correction factor by comparing the total (including central and satellite) halo mass functions from Behroozi et al. (2013) to those from the Bolshoi–Planck simulation (Rodríguez-Puebla et al. 2016), and find the best-fitting halo mass and redshift dependence of this correction. After this correction, the halo statistics used in this work are suitable for studying the evolution of haloes from $10^{10}$ to $10^{15} M_\odot$.

Haloes grow by smooth accretion and halo mergers. For average halo mass accretion histories, we use the fitting formulae in appendix H of Behroozi et al. (2013). For halo mergers, we fit merger rates from the UNIVERSEMACHINE (Behroozi et al. 2019), with full details and formulae in appendix B of Zhang et al. (2023c).

In addition to statistical halo assembly histories, we also use individual haloes from the MultiDark Planck 2 simulation (MDPL2; Klypin et al. 2016) to study the evolution of $z = 6$ haloes hosting SMBHs with $M_\bullet > 10^9 M_\odot$. The MDPL2 simulation also assumes Planck cosmology (Planck Collaboration XIII 2016), and has a simulation box size of 1 Gpc $h^{-1}$. The mass of each dark matter particle is $1.51 \times 10^9 M_\odot h^{-1}$, which enables the simulation to resolve haloes down to $M_{\text{peak}} \sim 10^{11} M_\odot$. As we will show in this work, the combination of the mass resolution and simulation box volume makes MDPL2 suitable for studying the host halo growth histories of $M_\bullet > 10^9 M_\odot$ SMBHs at $z \sim 6$, as well as the descendants of typical haloes hosting GNz-11 and CEERS_1019.

## 3.2 Observational data constraints

To constrain the 54 free parameters in TRINITY, we have originally compiled a diverse data set including galaxy and SMBH observables. The galaxy data include: galaxy stellar mass functions (SMFs), quenched fractions (QFs), cosmic star formation rates (CSFRs), average specific star formation rates (SSFRs), and UV luminosity functions (UVLFs). Collectively, these data cover a redshift range of $0 \le z \le 10$. The SMBH observables are: quasar luminosity functions (QLFs), quasar probability distribution functions (QPDFs), active black hole mass functions (ABHMFs), the $z = 0 M_\odot - M_{\text{bulge}}$ relation, and the observed $M_\bullet$ distribution of high redshift bright quasars. These SMBH observables collectively cover a redshift range of $0 \le z \le 6.5$. These datasets are summarized in tables 3–10 of Zhang et al. (2023c), and we refer readers to appendix C of Behroozi et al. (2019) and appendix D of Zhang et al. (2023c) for the details about the homogenization of galaxy and AGN data, respectively.

Since its launch, JWST has found more high-redshift UV-bright galaxies than predicted by theoretical models such as the UNIVERSEMACHINE. A higher number density of high-redshift galaxies likely implies more high-redshift SMBHs, which would systematically change our prediction for the early Universe. We thus added the $9 \lesssim z \lesssim 13$ galaxy UVLFs from Harikane et al. (2023b) to our data compilation, which includes only galaxies confirmed with spectroscopy. 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). For full details of calculating galaxy UV magnitudes, we refer the readers to appendix B of Zhang et al. (2023c).

Adding $9 \lesssim z \lesssim 13$ galaxy UVLFs does not qualitatively change in the results presented by Zhang et al. (2023c), Zhang et al. (2023a), or Zhang et al. (2023b). But this addition increases the number of UV-bright galaxies at $z \gtrsim 9$ and thus the galaxy masses at fixed halo masses, leading to more massive SMBHs in the same haloes. As shown in Fig. 1, this increase in galaxy number density is more significant at $z > 10$. This is because at $z \lesssim 10$, the new JWST UVLFs (solid circles) are largely consistent with existing UVLFs constraints (solid stars). But without $z > 10$ constraints, the original TRINITY underproduced UV-bright galaxies, 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 \lesssim z \lesssim 13$ galaxy UVLFs to our data constraints. We have also updated the plots in previous TRINITY papers (0:bold) here (0:bold).
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The comparison 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 (circles; Harikane et al. 2023b). The difference between the original and updated TRINITY models increases with redshift. See Section 3.2.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The observed median $M_{\bullet}-M_{\star}$ relations for bright quasars selected by different bolometric luminosity thresholds at $z=6-10$. The dashed, dash-dotted, dotted, and dash-double-dotted curves are the $M_{\bullet}-M_{\star}$ relations for bright quasars above different luminosity limits, which are annotated next to each line. These $M_{\bullet}-M_{\star}$ relations are nearly redshift-independent, so we only show the ones at $z=10$ for visual clarity. The shaded regions are the $1-\sigma$ lognormal scatter (~0.6 dex) around the $M_{\bullet}-M_{\star}$ relations for bright quasars, which includes the intrinsic scatter around the $M_{\bullet}-M_{\star}$ relation (~0.3 dex), as well as the typical random scatter in $M_{\bullet}$ from virial estimates (~0.5 dex). Note that, the shaded regions do not include the model uncertainties mentioned at the end of Section 2.2. Since the magnitude of this scatter is nearly luminosity-independent, we only show it for the brightest and faintest quasars for clarity. The solid lines of different colours represent the intrinsic median $M_{\bullet}-M_{\star}$ relations for all SMBHs at each redshift. UHZ1 (Bogdan et al. 2023; Goulding et al. 2023), CEERS.00717 (Harikane et al. 2023a), CEERS.1019 (Larson et al. 2023) and GN-z11 (Maiolino et al. 2023a) are shown in the diamond, the triangle, the pentagon, and the star, respectively. $M_{\bullet}$ uncertainties of 0.5 dex are assumed for CEERS.00717, CEERS.1019, and GN-z11. For UHZ1, the upper uncertainty of $M_{\bullet}$ is assumed to be 0.5 dex, but the lower uncertainty corresponds to changing the assumed Eddington ratio from 1 to 5, i.e. 0.7 dex. The galaxy mass uncertainties for these four high-redshift AGNs are taken from the source papers. We also show the $z \sim 6$ bright quasars with median log $L_{\text{bol}} \sim 47$ compiled by Izumi et al. (2021) in solid circles with error bars (~0.5 dex). All the observed data points are at least broadly consistent with TRINITY predictions, except for UHZ1. See Section 4.1.</figcaption></figure>

# 4 RESULTS

We present the offset between the $z=6-10 M_{\bullet}-M_{\star}$ relations for high redshift bright quasars vs. all SMBHs in Section 4.1, and the cosmic abundance of $M_{\bullet} > 10^{9} M_{\odot}$ and $10^{10} M_{\odot}$ SMBHs in Section 4.2. In Section 4.3, we show the $z=0$ descendant mass distribution of $z=6$ host haloes of $M_{\bullet} > 10^{9} M_{\odot}$ and $10^{10} M_{\odot}$ SMBHs, and discuss the prospect of detecting these descendants in the local Universe.

## 4.1 The biased $M_{\bullet}-M_{\star}$ relation for bright quasars between $6 \le z \le 11$

In Zhang et al. (2023a), we showed the median $M_{\bullet}-M_{\star}$ relation for bright quasars at $z \sim 6$, and concluded that overmassive SMBHs compared to the observed local $M_{\bullet}-M_{\star}$ relation are consistent with arising purely from Lauer bias against less massive AGNs (Lauer et al. 2007). In Fig. 2, we extend our predictions to even higher redshifts, up to $z=10$. The solid lines denote the intrinsic (i.e. including active and inactive SMBHs) median $M_{\bullet}-M_{\star}$ relation at different redshifts. The non-solid curves are the $M_{\bullet}-M_{\star}$ relations for quasars above various bolometric luminosity limits, which are labelled next to each line. Given their weak redshift dependence, we only show these quasar $M_{\bullet}-M_{\star}$ relations at $z=10$ for visual clarity. The shaded regions correspond to $1-\sigma$ lognormal scatter of ~0.6 dex around these quasar $M_{\bullet}-M_{\star}$ relations, which includes: (1) the intrinsic scatter around the $M_{\bullet}-M_{\star}$ relation (~0.3 dex, according to the best-fitting TRINITY model); and (2) the random scatter around $M_{\bullet}$ from virial estimates (0.5 dex; see e.g. McLure & Dunlop 2004; Vestergaard & Osmer 2009). The scatter is nearly luminosity-independent, so we only show the scatter for the brightest quasars with log $L_{\text{bol}} > 48$ and fainter AGNs with log $L_{\text{bol}} > 45$.

Between $6 \le z \le 10$, the $M_{\bullet}-M_{\star}$ relations for bright quasars continue to be higher than the intrinsic relation for all SMBHs. This is because at high redshifts, different SMBHs share very similar Eddington ratio distributions (see Section 2.1). As such, more massive SMBHs are brighter and over-represented in flux- or luminosity-limited samples. At fixed luminosity limits, the $M_{\bullet}-M_{\star}$ relations for quasars barely evolve with redshift, despite the evolution

in the intrinsic $M_{\bullet}-M_{\star}$ relation. This is because at fixed SMBH mass, the Eddington ratio distribution evolves very mildly at $z \gtrsim 6$. Thus, selecting quasars by their luminosities yields similar SMBH masses across different redshifts.

We also note that the typical $M_{\bullet}$ for quasars is essentially independent of host galaxy mass at log $M_{\star} \lesssim 9.5$. This is because high-redshift SMBHs grow at average Eddington ratios slightly below the Eddington rate, and hyper-Eddington AGNs are very rare. This upper limit in Eddington ratios leads to a lower limit in SMBH mass for luminosity-limited samples. For the brightest quasars (e.g. log $L_{\text{bol}} \ge 48$), the typical SMBH mass remains above $3 \times 10^8 M_{\odot}$ even for galaxies below the same mass. Hence, the detection of any log $L_{\text{bol}} \ge 48$ quasars in a $M_{\star} < 3 \times 10^8 M_{\odot}$ galaxy implies: (1) it has an Eddington ratio much higher than predicted by TRINITY; and/or (2) it is in the outsized black hole galaxy (OBG) phase, where the SMBH could be more massive than the galaxy (Agarwal et al. 2013; Natarajan et al. 2017).

In this work, we assume that every AGN above those luminosity limits will be detected and measured accurately. In other words, we ignore any potential biases due to the light from host galaxies. This assumption is more likely to hold for brighter quasars. For fainter AGNs, e.g. log $L_{\text{bol}} \sim 45$, contamination from host galaxies
makes it more difficult to measure SMBH masses, even if the AGN luminosities are by themselves above the limits. As a result, overmassive SMBHs would still be preferred (see, e.g. Volonteri, Habouzit & Colpi 2023). It takes detailed modeling of high-redshift galaxy and SMBH spectral energy distributions (SEDs) to properly account for this bias, which we defer to a future work.

Besides Lauer bias, Li et al. (2024) showed that galaxy mass measurement uncertainties and the selection based on broad emission line widths would further bias the $M_{\bullet}-M_{*}$ relation. Therefore, our predicted $M_{\bullet}-M_{*}$ relations for fainter AGNs (e.g. with $\log L_{\text{bol}} > 45$) should be regarded as lower limits.

We also show the four new highest-redshift AGNs detected by JWST, CEERS.00717 (Harikane et al. 2023a; the triangle), CEERS.1019 (Larson et al. 2023; the pentagon), UHZ1 (Bogdan et al. 2023; Goulding et al. 2023; the diamond) and GN-z11 (Maiolino et al. 2023a; the star) in Fig. 2. Since the SMBH masses of CEERS.00717, CEERS.1019, and GN-z11 are estimated with virial estimates, we adopt an $M_{\bullet}$ uncertainty of 0.5 dex for these three AGNs. The $M_{\bullet}$ of UHZ1 was estimated from its X-ray luminosity assuming Eddington-limited accretion. Given that SMBHs accreting at five times the Eddington rate have been detected at similar redshifts (GNz-11), we adopt a $M_{\bullet}$ uncertainty of 0.7 dex (i.e. a factor of 5) to reflect the uncertainty in Eddington ratio. Galaxy mass uncertainties are directly taken from these source papers. In several cases, stellar masses have been determined without including an AGN component in the fit, and the uncertainties in the star formation histories can induce large systematic uncertainties (see e.g. the summary in table 1 of Volonteri et al. 2023).

Compared to the brighter quasars with median $\log L_{\text{bol}} \sim 47$ compiled by Izumi et al. (2021), all these JWST AGNs have lower SMBH masses. This is consistent with TRINITY's prediction of Lauer bias in the $M_{\bullet}-M_{*}$ relation. At face value, the SMBH in CEERS.00717 seems much more massive than TRINITY's prediction for faint AGNs. But given the substantial uncertainty in galaxy mass, we decide to defer conclusions on CEERS.00717 until more precise measurements of galaxy mass are available. The mass and luminosity measurements of CEERS.1019 and GN-z11 are broadly consistent with TRINITY predictions, lying above the relation for $\log L_{\text{bol}} > 45$ at $\sim 1.5\sigma$ and $\sim 0.8\sigma$, respectively. This could be due to: (1) small number statistics; (2) selection effects beyond Lauer bias (Volonteri et al. 2023; Li et al. 2024); and/or 3) TRINITY underpredicting SMBH masses compared to the real Universe. Further observations and better characterizations of all selection biases are needed to ascertain the nature of this difference. In Section 5.2, we compare CEERS1019 and GN-z11 with other predictions from TRINITY, such as AGN Eddington ratios and number densities.

Finally, we make a few remarks about UHZ1. Assuming Eddington accretion, the $M_{\bullet}$ of UHZ1 is more than 1 dex higher than expected for AGNs above similar luminosities ($\log L_{\text{bol}} \ge 45.7$, the dash-dotted line). But another super-Eddington AGN has been found at $z \gtrsim 10$, i.e. GN-z11. If UHZ1 is actually also accreting at a similar Eddington ratio, the resulting $M_{\bullet}$ would be marginally consistent with TRINITY predictions. A caveat of this argument is that the current estimate of UHZ1's luminosity (and $M_{\bullet}$) is only a $1\sigma$ lower limit from X-ray spectral fitting. This value was reported instead of the actual best-fitting value due to the strong degeneracy between the intrinsic X-ray luminosity and gas column density used in the spectral fitting process. Thus, choosing higher intrinsic luminosity values will increase the tension between UHZ1 and TRINITY predictions.

The detection of UHZ1 itself also indicates a strong deviation from TRINITY predictions, if $M_{\bullet}$, $M_{*}$, and AGN luminosity measurements are taken at face value. TRINITY predicts that the probability of

detecting at least one UHZ1-like object<fnref n="1" /> in the UNCOVER field is $\sim 2 \times 10^{-7}$ (see equations 4 and 5). In this prediction, we assume $\sigma_{M_{\bullet}} = \sigma_{M_{*}} = 0.5$ dex rather than the values from Bogdan et al. (2023) and Goulding et al. (2023), to account for measurement uncertainties from SED fitting and virial estimates. Such a small probability results from: (1) the overmassiveness of the SMBH compared to the host galaxy, and (2) the extremely low-AGN duty cycles for low-mass galaxies (see Section 5.2.3). Hence, if the measurements of UHZ1 are accurate, its detection would be in stark contrast with TRINITY predictions. The large $M_{\bullet}/M_{*}$ ratio implies that UHZ1 may be experiencing an OBG phase, which is not accounted for in the current TRINITY model. If confirmed, AGNs like UHZ1 will provide critical information for future versions of TRINITY to account for this OBG scenario. One possible solution is to introduce a more sophisticated feedback scenario at high redshift that gives priority to SMBH formation followed by a brief phase of boosted star formation via SMBH-induced positive feedback in the emerging galaxy (Silk et al. 2024).

## 4.2 $M_{\bullet} > 10^9 M_{\odot}$ and $10^{10} M_{\odot}$ SMBHs at $z > 6$

The top panel of Fig. 3 shows the redshift density, $dN/dz$, of $M_{\bullet} > 10^9 M_{\odot}$ and $10^{10} M_{\odot}$ SMBHs in the whole observable Universe. The solid lines represent all (i.e. both active and dormant) SMBHs with intrinsic masses above these thresholds. From $z = 10$ to $z \sim 1–2$, $dN/dz$ increases rapidly with time due to high accretion rates. Below $z \sim 1–2$, $dN/dz$ decreases due to the decline of both SMBH growth and observable volume. The middle panel of Fig. 3 shows cumulative numbers, $N(z)$, where $N(z) = \int_{z}^{\infty} dN/dz dz$, of $M_{\bullet} > 10^9 M_{\odot}$ and $10^{10} M_{\odot}$ SMBHs. Beyond the local Universe, virial estimates remain the only way to measure black hole masses for large samples. These methods result in random scatter of 0.4–0.6 dex in observed black hole masses around the intrinsic values (McLure & Jarvis 2002; Vestergaard & Peterson 2006). This scatter causes an Eddington bias (Eddington 1913), which elevates the observed number of massive black holes above the true number (dashed curves in Fig. 3, Top and Middle Panels). Overall, an observational scatter of 0.5 dex elevates observed massive black hole number by a factor of $\sim 3–40$. This boost in black hole numbers is stronger for more massive objects and/or lower redshifts. This is because the massive end of the SMBH mass function becomes steeper towards lower redshifts (Fig. 4). The steeper mass function leads to a larger population of upscattered $M_{\bullet}$ measurements, and thus a stronger boost in the numbers of large SMBHs. With this Eddington bias, the cumulative numbers of $M_{\bullet} > 10^9$ and $10^{10} M_{\odot}$ SMBHs increase by five orders of magnitude from $z \sim 10$ to $z \sim 2$.

According to TRINITY, a $10^9 M_{\odot}$ SMBH would be $\sim 3–4$ dex above the median $M_{\bullet}-M_{*}$ relation at $z = 10$ (Fig. 2). The rarity of such SMBHs is also shown in the middle panel of Fig. 3: we expect to see only $\lesssim 100$ such SMBHs at $z \gtrsim 10$ over the whole sky, whether they are active or not. Such rarity results from the extreme difficulty of growing a $10^9 M_{\odot}$ SMBH by $z \sim 10$, which requires an appropriate combination of SMBH seed mass, Eddington ratios, duty cycle, and AGN efficiency (Pacucci & Loeb 2022). Therefore, the detection of massive BHs at $z \gtrsim 10$ will strongly constrain SMBH seeding mechanisms and their early growth histories. We will further compare the possible progenitor masses and growth speeds of $10^9 M_{\odot}$ SMBHs from Pacucci & Loeb (2022) and TRINITY in Section 5.3.

<footnote n="1">¹i.e. $\log M_{\bullet} \ge 7.5$ SMBHs brighter than $\log L_{\text{bol}} \ge 45.7$ and hosted by galaxies with $\log M_{*}, \text{obs} \sim 8.15$.</footnote>
<figure><img src="image_4.png" /></figure>

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

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

**Figure 3.** Top Panel: The number of very massive black holes per unit redshift in the whole observable universe (i.e. dN/dz). The solid lines show results assuming no observational scatter in the measured black hole masses, and the dashed lines show results assuming an observational scatter of σvir = 0.5 dex, typical for virial estimates (see Section 4.2 for details). The dotted lines show the numbers of SMBHs based on the luminosity-dependent selection function shown in the bottom panel. Middle Panel: The cumulative number of very massive black holes as a function of redshift. The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. These intervals reflect the ranges that the model predictions are allowed to vary, while keeping the underlying SMBH properties matching input data constraints within error bars (also see Section 2.2). Bottom Panel: AGN completeness, i.e. the fraction of quasars that would survive the combination of dropout and colour selection, as a function of bolometric luminosity and redshift from Tee et al. (2023) and Euclid Collaboration (2019). The number of SMBHs with M$_{\bullet} > 10^9 M_{\odot}$ in the observable universe increases by five orders of magnitude from z ~ 10 to z ~ 2, and by another factor of ~3 from z ~ 2 to z = 0 See Section 4.2.

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

**Figure 4.** Total SMBH mass functions (BHMFs) from $z = 0–15$. The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. These intervals reflect the ranges that the model predictions are allowed to vary, while keeping the underlying SMBH properties matching input data constraints within error bars (also see Section 2.2). The BHMF decreases strongly in normalization and flattens towards higher redshifts. See Section 4.2.

To better understand the SMBH detection capability of future wide field surveys, we also made predictions of expected SMBH numbers from the *Legacy Survey of Space and Time (LSST)* and the *Euclid* wide survey. To calculate the number of detectable quasars at different luminosities and redshifts, we calculate the AGN completeness, i.e. the fraction of quasars that would survive the combination of dropout and color selections. First, we construct a uniform magnitude-redshift grid for simulating quasar photometry. Secondly, we use quasar observed continuum slope and emission line strengths to build empirical quasar spectra, with the assumption that quasar SEDs are redshift independent (Fan 1999). The distribution of continuum slopes, emission line features (full widths at half maximum and equivalent widths), and quasar spectral characteristics (Baldwin effect and blueshifted lines) are designed to match the SDSS BOSS quasars in the redshift range $2.2 < z < 3.5$ (Ross et al. 2012). This framework has been extended to include quasars at higher redshifts ($z = 4–10$) in a number of recent works for studying quasar selection functions (McGreer et al. 2013; Jiang et al. 2016; Yang et al. 2016; McGreer et al. 2018; Wang et al. 2019; Tee et al. 2023). The quasar selection function is a joint colour selection modified from Richards et al. (2002) and Euclid Collaboration (2019) at $z = 0–6$ and $z = 6–10$, respectively, to maximize the detection of SMBHs. In the low-$z$ regime, the original selection is not sensitive to faint quasars and results in a deficit from $z = 0–3$ at the faint end. To recover faint populations, we add an additional colour selection of $u - g > 0.5$ and $g - r > 0.1$, missed by the original selection. The modified selection completeness is similar to Richards et al. (2002) but reaches $> 3$ mag deeper. In the high-$z$ regime, we are able to select quasars with a combination of dropout and colour cuts, i.e. $z - Y > 1.5$ for $z = 6–9$ and $Y - J > 1.5$ for $z = 9–10$. Here, the selection completeness is slightly deeper compared with Euclid Collaboration (2019), as the latter focused on minimizing contaminants. In the bottom panel of Fig. 3, we show quasar completeness as a function of bolometric luminosity and redshift, and the differential (cumulative) number of detectable massive black holes is shown in dotted lines in the top (middle) panel of Fig. 3.

At $z \gtrsim 4$, the vast majority of SMBHs with $M_{\bullet, \text{obs}} > 10^9 M_{\odot}$ are active, so the bulk of these SMBHs are above the flux limits of LSST and Euclid. We note that the number of SMBHs in Fig. 3 include
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Top Panel: The host halo mass functions of SMBHs with intrinsic (solid lines) and observed (dashed lines) $M_{\bullet} > 10^9 M_{\odot}$ and $M_{\bullet} > 10^{10} M_{\odot}$ at $z=6$, calculated by TRINITY. We also show the $z=6$ halo mass function from the MDPL2 simulation in solid line, scaled down to 10 per cent of the original value for ease of comparison. Bottom Panel: The $z=0$ descendant halo mass functions of the $z=6$ haloes hosting SMBHs with intrinsic (solid lines) and observed (dashed lines) $M_{\bullet} > 10^9 M_{\odot}$ and $M_{\bullet} > 10^{10} M_{\odot}$. They are calculated by following individual haloes in the MDPL2 simulation across time, and weighting each of them by the ratios between the host halo mass functions and the total halo mass functions in the top panel. We also show the $z=0$ halo mass function in the MDPL2 simulation in solid line, scaled down by a factor of 100. $M_{\bullet} > 10^9 M_{\odot}$ at $z=6$ are hosted by ~$10^{12}-10^{13} M_{\odot}$ haloes, which can grow into $10^{14}-10^{15} M_{\odot}$ by $z=0$. But at $z=0$, they only make up ~1 per cent of all the massive haloes. See Section 4.3.</figcaption></figure>

both obscured and unobscured SMBHs. Given the uncertainties in the obscured AGN fraction towards higher redshifts, we opt not to show results for unobscured AGNs in Fig. 3, but only give a rough estimate: at $L_{bol} \gtrsim 10^{46}$ erg s⁻¹, we expect ~40 per cent (or 80 per cent) of such SMBHs to be unobscured, based on the obscured fractions as functions of AGN X-ray luminosity given by Merloni et al. (2014) (or Ueda et al. 2014). Below $z \sim 4$, massive black holes experience a significant decline in activity level. Consequently, fewer and fewer massive black holes at lower redshifts will be detectable with LSST.

## 4.3 Host halos of massive black holes at $z=6$ and their descendants at $z=0$

We show the host halo mass functions (HMFs) of $M_{\bullet} > 10^9$ and $10^{10} M_{\odot}$ SMBHs at $z=6$, in the top panel of Fig. 5 (solid curves). When using virial estimates to measure $M_{\bullet}$, there is a random scatter

in the measured $M_{\bullet}$ around the intrinsic values. As such, we also show the host halo mass functions for observed $M_{\bullet} > 10^9$ and $10^{10} M_{\odot}$ in dashed curves as well. In this comparison, we assume a random scatter of 0.5 dex in measured $M_{\bullet}$. This random scatter induces Eddington bias, which enhances the host HMF normalizations and shifts HMFs towards lower masses. The median host halo mass for SMBHs with intrinsic $M_{\bullet} > 10^9 M_{\odot}$ is ~3 × $10^{12} M_{\odot}$, compared to ~2 × $10^{12} M_{\odot}$ for SMBHs with observed $M_{\bullet} > 10^9 M_{\odot}$. Such halo masses are consistent with the findings of Shimasaku & Izumi (2019), where halo masses were measured with the width of [C II] 158 µm lines. For $M_{\bullet} > 10^{10} M_{\odot}$ SMBHs, their host haloes are typically ~5 × $10^{12} M_{\odot}$ and ~3 × $10^{12} M_{\odot}$ for intrinsic and observed $M_{\bullet}$, respectively. At $z=6$, host haloes of $M_{\bullet} > 10^9$ and $10^{10} M_{\odot}$ SMBHs only make up ~10 per cent and 1 per cent of all the haloes at their respective median host halo masses, which correspond to ~4.5–5.5σ peaks in Gaussian random fields (Eisenstein & Hu 1998; Diemer 2018). In such haloes, SMBH could grow via cold flows (e.g. Feng et al. 2014 and Latif et al. 2022), accretion triggered by protogalaxy mergers (e.g. Li et al. 2007), or a combination of the two (e.g. Dubois et al. 2012). Observationally, we expect to see overdensities of galaxies around the SMBHs in these massive haloes (Wang et al. 2023). In future work, we will create realistic mock galaxy + AGN catalogues to quantify the typical magnitude and stochastic variance of galaxy overdensities around $z \gtrsim 6$ quasars, and compare them with results from, e.g. the ASPIRE (Wang et al. 2023) and EIGER (Kashino et al. 2023) surveys.

Ever since the discovery of SMBHs with $M_{\bullet} > 10^9 M_{\odot}$, there have been discussions about the potential $z=0$ descendants of these early SMBHs (e.g. Angulo et al. 2012; Fanidakis et al. 2013; Tenneti et al. 2018). To address this question with TRINITY, we follow each $z=6$ halo in the MDPL2 simulation down to $z=0$, and examine their descendant halo masses. We weight each progenitor- descendant pair by the fraction of $z=6$ haloes that would host a $>10^9$ and $10^{10} M_{\odot}$ SMBH as a function of progenitor mass (i.e. the ratios between the curves for $>10^9$ and $10^{10} M_{\odot}$ SMBHs and the total HMF in Fig. 5, Top Panel). As shown in the bottom panel of Fig. 5, the descendant HMFs are much broader than the progenitor HMFs due to the diversity in individual halo assembly histories. At $z=0$, descendants with $M_{\bullet}(z=6) > 10^9 M_{\odot}$ typically have $M_{peak} \sim 3-7 \times 10^{14} M_{\odot}$. However, these haloes only make up ~1 per cent of the haloes at similar mass scales (the orange solid line). Given the small comoving volume in the local Universe, there is only ~1 such relic SMBH closer than $z=0.05$. Local SMBHs with $M_{\bullet}(z=6) > 10^{10} M_{\odot}$ have similar host halo masses, but they are much rarer, by a factor of ~30–100.

## 5 DISCUSSION

In Section 5.1, we discuss the physical implications of the steeper intrinsic $M_{\bullet}-M_{*}$ relation between $7 \lesssim z \lesssim 10$ from TRINITY. In Section 5.2, we compare CEERS.1019 and GN-z11 with more TRINITY predictions at $6 \lesssim z \lesssim 11$. In Section 5.3, we compare the progenitor masses and growth rates of $\gtrsim 10^9 M_{\odot}$ SMBHs at $z \sim 9$ from Pacucci & Loeb (2022) with those from TRINITY. Section 5.4 discusses the constraining power of $z=6-10$ quasar luminosity functions on the halo–galaxy–SMBH connection.

### 5.1 Physical Implications of the $M_{\bullet}-M_{*}$ relation at $z \gtrsim 6$ from TRINITY

According to TRINITY, the intrinsic $M_{\bullet}-M_{*}$ relation between $7 \lesssim z \lesssim 10$ is mildly evolving, and the relation is steeper towards higher
redshifts (see also Volonteri & Stark 2011). This implies that SMBH growth is slower than host galaxy growth in $M_* < 10^{10}$ M$_\odot$ galaxies in the early Universe. Several physical mechanisms could contribute to such limited growth: (1) strong supernova feedback in low-mass galaxies preventing efficient SMBH accretion (as in, e.g. Dubois et al. 2015; Bower et al. 2017; Pillepich et al. 2018); (2) erratic motions of SMBHs due to shallow potential wells (as in, e.g. Habouzit, Volonteri & Dubois 2017; Bellovary et al. 2019; Dekel, Lapiner & Dubois 2019, though c.f. Choksi et al. 2017); and (3) Feedback from SMBHs themselves (e.g. as in ROMULUS; Tremmel et al. 2017). Limited SMBH growth at $z \sim 6-7$ also implies that AGN are unlikely to dominate the cosmic reionization of intergalactic medium, as also suggested by e.g. Finkelstein et al. (2019).

## 5.2 CEERS_1019 and GN-z11 in the context of TRINITY

### 5.2.1 SMBH mass and Eddington ratios

Fig. 6 shows the $z = 10$ bolometric quasar luminosity function (QLF) decomposed into bins of SMBH mass ($M_\bullet$, top panel) and Eddington ratio ($\eta_{\text{rad}}$, bottom panel), predicted by TRINITY. At $z \sim 10$, the QLF peaks at around $10^{45}$ erg s$^{-1}$. In addition, SMBHs with $6 < \log M_\bullet < 7$ (top panel) and $0 < \log \eta_{\text{rad}} < 1$ contribute significantly the QLF at around $10^{45}$ erg s$^{-1}$. Given that both GN-z11 and CEERS_1019 have bolometric luminosities of $\sim 10^{45}$ erg s$^{-1}$, their masses and Eddington ratios (solid and dashed vertical lines) are consistent with TRINITY's predictions. We also note that in TRINITY, $z \sim 10$ SMBHs grow at slightly below the Eddington rate, when averaged over both their active and dormant phases. Consequently, although GN-z11 is estimated to be accreting at a super-Eddington rate, such fast accretion phases are likely intermittent. Two objects are far from a statistically representative sample, so more detections are needed to further confirm/disprove TRINITY's predictions of typical SMBH mass and Eddington ratios for $z \gtrsim 9$ AGNs.

### 5.2.2 Abundance of $M_\bullet > 10^6/10^7$ M$_\odot$ SMBHs between $6 < z < 11$

Fig. 7 shows the cumulative number of SMBHs above $10^6$ M$_\odot$ and $10^7$ M$_\odot$ from $z = 6-15$. The solid and dashed curves represent the number of all (i.e. active + dormant) SMBHs with intrinsic and observed $M_\bullet$ above these mass limits, respectively. The dotted curves are similar to the dashed curves, but for SMBHs above the bolometric luminosity limit of $10^{45}$ erg s$^{-1}$. This limit has been chosen to match the luminosities of CEERS1019 and GN-z11. Similar to Fig. 3, these curves represent the expected SMBH numbers over the whole sky. Compared to more massive SMBHs in Fig. 3, Eddington bias causes a smaller boost in the number of SMBHs, i.e. a factor of only $\sim 1.5$. This is because the total BHMF is flatter at the low-mass end, which reduces the magnitude of Eddington bias. According to TRINITY, we expect to find $\sim 30000$ and $\sim 3000$ SMBHs with observed $M_\bullet > 10^6$ and $10^7$ M$_\odot$ at $z > 10$, respectively. However, only a small fraction of these SMBHs are bright enough to be detected. Given the nearly mass-independent Eddington ratio distributions from TRINITY, fewer low-mass SMBHs will survive the same luminosity limit. Consequently, at $L_{\text{bol}} > 10^{45}$ erg s$^{-1}$, the numbers of SMBH with $M_\bullet > 10^6$ and $10^7$ M$_\odot$ at $z > 9$ decrease to $\sim 1000$. The small number of $L_{\text{bol}} > 10^{45}$ erg s$^{-1}$ AGNs at $z > 9$ indicates a small probability of detecting GN-z11 and CEERS_1019 analogues in deep surveys, which we will further quantify in Section 5.2.3.

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

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Top Panel: The $z = 10$ bolometric quasar luminosity function (QLF) decomposed into bins of different SMBH mass ($M_\bullet$). Bottom Panel: The $z = 10$ bolometric quasar luminosity function (QLF) decomposed into bins of different SMBH Eddington ratio ($\eta_{\text{rad}}$). The solid and dashed vertical lines represent the bolometric luminosities of CEERS_1019 and GN-z11, respectively. The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. These intervals reflect the ranges that the model predictions are allowed to vary, while keeping the underlying SMBH properties matching input data constraints within error bars (also see Section 2.2). The SMBH masses and Eddington ratios of CEERS_1019 and GN-z11 are representative of the $L_{\text{bol}} \sim 10^{45}$ erg s$^{-1}$ AGNs in TRINITY. See Section 5.2.1.</figcaption></figure>

### 5.2.3 AGN fraction among $6 < z < 11$ galaxies

Fig. 8 shows the fraction of galaxies hosting AGNs with bolometric luminosities above $10^{45}$ erg s$^{-1}$, $f_{\text{gal}}$, as a function of galaxy mass and redshift. These fractions include both obscured (including Compton-thick obscured) and non-obscured AGNs. The stochastic uncertainties in $f_{\text{gal}}$ at $z = 6$ and $z = 11$ are shown in shaded regions. We also show the galaxy masses of GN-z11 and CEERS_1019 in red solid and dashed lines, respectively.

At $z \gtrsim 6$, $f_{\text{gal}}$ is a strong function of galaxy mass, increasing from $\lesssim 10^{-6}$ at $M_* \sim 10^9$ M$_\odot$ to $\sim 1$ at $M_* \sim 3 \times 10^{10}$ M$_\odot$. This is an extrapolation based on quasar probability distribution functions (QPDFs) from Aird, Coil & Georgakakis (2018), which requires very low-AGN duty cycles in low-mass galaxies. In addition, lower-mass galaxies host smaller SMBHs, which further decreases AGN incidence at the low-mass end. The high-AGN incidence in massive galaxies suggests that it is essential to properly account for any AGN contribution when trying to characterize high-redshift massive
<figure><img src="image_12.png" /><figcaption>**Figure 7.** The cumulative number of SMBHs with $M_{\bullet} > 10^6 M_{\odot}$ and $M_{\bullet} > 10^7 M_{\odot}$. The line styles are the same as in Fig. 3, but we adopt a constant luminosity limit of $10^{45}$ erg s$^{-1}$ in this figure. The darker shaded region is the 68th percentile range (from MCMC) for $M_{\bullet} > 10^6 M_{\odot}$ with no scatter in $M_{\bullet}$$_{,obs}$. The relative uncertainty in the number of $M_{\bullet} > 10^7 M_{\odot}$ SMBHs (red solid curve) is very similar to $M_{\bullet} > 10^6 M_{\odot}$ SMBHs, so we do not show its confidence intervals for visual clarity. and the lighter blue (red) shaded region is the the 68th percentile range for $M_{\bullet} > 10^6 M_{\odot} (10^7 M_{\odot})$ above the bolometric luminosity of $10^{45}$ erg s$^{-1}$, with $M_{\bullet}$$_{,obs} = 0.5$ dex. All the shaded regions reflect the ranges that the model predictions are allowed to vary, while keeping the underlying SMBH properties matching input data constraints within error bars (also see Section 2.2). The number of $M_{\bullet} > 10^6 M_{\odot}$ SMBHs increases by ~6 order of magnitudes, but only a small fraction of them are brighter than $10^{45}$ erg s$^{-1}$. See Section 5.2.2.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 8.** The fraction of galaxies that host AGNs above the bolometric luminosity of $10^{45}$ erg s$^{-1}$ as a function of galaxy mass ($M_{\star}$) and redshift. The solid and dashed vertical lines denote the estimated galaxy masses for Maiolino et al. (2023a) and Larson et al. (2023), respectively. For visual clarity, we only show the 68th percentile ranges (from MCMC) for $z = 6$ and $z = 11$ in shaded regions. These shaded regions reflect the ranges that the model predictions are allowed to vary, while keeping the underlying SMBH properties matching input data constraints within error bars (also see Section 2.2). $f_{\text{gal}}$ strongly increases towards higher galaxy mass and lower redshift. See Section 5.2.3.</figcaption></figure>

galaxies (Volonteri et al. 2017). On the flip side, it also means that detections of high-redshift massive galaxies that have not included AGN contributions likely need to be revisited (e.g. Labbé et al. 2023). In the low-mass regime, it is statistically more acceptable to ignore potential AGN contributions when characterizing galaxies.

However, it is still essential to account for AGN contributions when characterizing host galaxies of confirmed AGNs like GN-z11 and CEERS-1019.

The small AGN fraction among low-mass galaxies also implies that GN-z11 and CEERS_1019 are rare cases where low-mass galaxies do host bright AGNs, which seem to be in tension with observations given the survey areas of JADES and CEERS (see also Maiolino et al. 2023b for an estimate of AGN fraction). Based on equation 5, TRINITY predicts probabilities of ~0.04 and ~0.02 to detect at least one GN-z11 and CEERS_1019 analogue in GOODS-S + N and CEERS fields, respectively. For this calculation, we adopt $\sigma_{M_{\star}} = \sigma_{M_{\bullet}} = 0.5$ dex, as we did for UHZ1 in Section 4.1. Although these probabilities are still relatively small, they do show that the GN-z11 and CEERS_1019 are less unexpected than UHZ1. As we discuss in Section 5.2.4, the constraint from GN-z11 and CEERS_1019 is not strong enough to significantly change TRINITY model results. Hence, in addition to a better understanding of galaxy mass measurements and various selection effects, non-detection over a much bigger survey area or a larger low-mass galaxy sample with similar AGNs is needed to confirm or disprove high-redshift predictions from TRINITY.

### 5.2.4 Using CEERS_1019 and GN-z11 as constraints on TRINITY

In this section, we quantify the constraint by CEERS_1019 and GN-z11 on the current TRINITY model if they were added to the current data compilation. We do not include UHZ1 in this exercise because the only evidence for UHZ1's AGN nature is the X-ray detection, for which there are still mixed opinions on whether the X-ray photons come from a foreground source by coincidence. Given the uncertainties in their host galaxy masses, we chose not to compare them with the number density of AGNs hosted by similar-mass galaxies, as in Fig. 8. Instead, we converted the detection of CEERS_1019 and GN-z11 into the cumulative number density of AGNs above $L_{\text{bol}} = 10^{45}$ erg s$^{-1}$ between $z = 8-12$. To calculate the comoving volume, we combined the survey areas of CEERS, GOODS-S, and GOODS-N fields, which amount to 260 arcmin$^2$. The fiducial TRINITY model predicts ~0.15 AGNs above $L_{\text{bol}} = 10^{45}$ erg s$^{-1}$ between $z = 8-12$ in these fields. Using the two detected AGNs as a lower limit and assuming Poisson statistics, this corresponds to a $\chi^2_{z>9} \sim 9.2$. We added this additional number density constraint to TRINITY, and reran the MCMC. We found that $\chi^2_{z>9}$ is too small to significantly change the posterior distribution of model parameters, compared to the total $\chi^2_{\text{tot}} \sim 770$ from other data. In other words, this constraint in number density is not strong enough to change the empirical galaxy-SMBH connection that is already constrained by all the other existing galaxy and SMBH data. This further confirms the need for either non-detections of such high-redshift AGNs over a larger volume and better host galaxy measurements to confirm TRINITY's predictions, or more detections per unit volume to provide strong enough constraints on TRINITY.

If future observations (with selection effects corrected) confirmed that the AGN fraction is indeed much higher for $M_{\star} \lesssim 10^9 M_{\odot}$ galaxies than predicted by TRINITY (Fig. 8), then TRINITY would likely have underestimated the following SMBH properties: (1) the typical SMBH masses; (2) typical AGN Eddington ratios; and/or (3) AGN duty cycles in these low-mass galaxies. Physically, this implies that high-redshift SMBHs in low-mass galaxies likely grow faster and more consistently than their low-redshift counterparts, which is not captured by extrapolating low-redshift data with the
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 9.** Top Panel: The host halo mass functions of SMBHs with intrinsic (solid curves) and observed (dashed curves) $M_{\bullet} > 10^6 M_{\odot}$ at $z = 10.6$ and $M_{\bullet} > 10^7 M_{\odot}$ at $z = 8.7$. Bottom Panel: The $z = 0$ descendant halo mass functions of halos hosting intrinsic (solid curves) and observed (dashed curves) $M_{\bullet} > 10^6 M_{\odot}$ at $z = 10.6$ and $M_{\bullet} > 10^7 M_{\odot}$ at $z = 8.7$. We also show the total host halo mass functions $z = 0$ in solid curve, scaled down by a factor of 100. $M_{\bullet} > 10^6 M_{\odot}$ at $z > 9$ are hosted by $\sim 10^{11}-10^{12} M_{\odot}$ haloes, which can grow into $10^{13}-10^{15} M_{\odot}$ by $z = 0$. But at $z = 0$, they only make up $\sim 1$ per cent of all the massive haloes. See Section 5.2.5.</figcaption></figure>

current parametrization. One possible solution to this tension is to introduce new model parameters to describe the high-redshift behaviour of these SMBH properties that are not captured by the current extrapolation. We will defer this exploration until we have a sufficient number of $z \gtrsim 8$ AGNs with robust measurements of host galaxy properties.

### 5.2.5 Host haloes of CEERS_1019 and GN-z11

In the top panel of Fig. 9, we show the host halo mass functions of SMBHs with intrinsic (solid curves) and observed (dashed curves) $M_{\bullet} > 10^6 M_{\odot}$ at $z = 10.6$ and $M_{\bullet} > 10^7 M_{\odot}$ at $z = 8.7$. These redshifts and mass limits are chosen to roughly match the observed properties of CEERS_1019 and GN-z11. Similar to the host halo mass functions of $M_{\bullet} > 10^9 M_{\odot}$ at $z = 6$ (Section 4.3, Fig. 5), the difference between the solid and dashed curves comes from the Eddington bias. GN-z11 is likely hosted in a $(2-3) \times 10^{11} M_{\odot}$ halo, whereas CEERS1019's typical host halo mass is $(3-5) \times 10^{11} M_{\odot}$. These host halo masses correspond to $\gtrsim 7\sigma$ peaks in Gaussian random fields at these redshifts, making them candidates for protocluster cores. This also means that the typical host haloes for CEERS_1019 and GN-z11 are rarer than

the luminous quasars at $z \sim 6$ at their redshifts. Our estimated host halo mass for GN-z11 is $\sim 1$ dex higher than the one by Scholtz et al. (2023). This discrepancy is likely due to the difference in how host halo masses are estimated: Scholtz et al. made this estimate based on the stellar mass measurement by Tacchella et al. (2023) and the stellar mass-halo mass relation from Behroozi et al. (2019), assuming halo mass-independent AGN duty cycles. In contrast, TRINITY predicts that AGN duty cycles are lower in low-mass haloes. In other words, low mass haloes are much less likely to host AGNs like GN-z11. These predictions are driven by extrapolating the QPDFs from Aird et al. (2018) to high redshifts. Therefore, TRINITY finds more massive host haloes for GN-z11, despite the very similar galaxy mass-halo mass relations at $z \gtrsim 10$.

The bottom panel of Fig. 9 shows the $z = 0$ descendant halo mass functions of haloes hosting $M_{\bullet} > 10^6 M_{\odot}$ at $z = 10.6$ and $M_{\bullet} > 10^7 M_{\odot}$ at $z = 8.7$. These descendant halo mass functions are calculated using the method introduced in Section 4.3. The typical descendant halo masses are $\sim (1-3) \times 10^{14} M_{\odot}$, but the distribution is broad due to the diversity in individual halo assembly histories. The descendants make up only $\lesssim 1$ per cent of all the haloes at similar mass scales. This makes them too rare to find in the local Universe. Our descendant mass scale is lower by $\sim 1$ dex compared to the descendant halo mass of GN-z11 obtained by Scholtz et al. (2023), which was based on the Millenium Simulation (Springel et al. 2005; Chiang, Overzier & Gebhardt 2013). This is because Scholtz et al. (2023) chose to look at the average $z = 10.6$ progenitor halo masses of $z = 0$ haloes above $10^{15} M_{\odot}$, and found that the average progenitor halo mass coincides with the halo mass estimate of GNz-11. Due to the scatter in individual halo assembly histories, this is different from looking at the $z = 0$ descendant halo masses of GNz-11 host haloes, which will lead to a smaller $z = 0$ halo mass. In this work, we chose the latter method because it is more representative of the expected descendant masses; the Scholtz et al. (2023) method is more appropriate for finding $z = 0$ objects whose typical progenitors look like GNz-11.

### 5.2.6 Potential descendants and progenitors of CEERS1019 and GN-z11

To further understand the potential progenitors and descendants of GN-z11 and CEERS_1019, we calculated the average$^2$ Eddington ratio as a function of $M_{\bullet}$ and time, and the average mass growth curves (Fig. 10) for GN-z11, CEERS_1019, and a seed of 100 $M_{\odot}$ at $z \sim 15$ (solid curves and circles). The shaded regions represent the statistical uncertainties from the MCMC, which reflects the uncertainty in the input observational constraints (see Section 2). In other words, they quantify the uncertainty of the average mass growth curves given all existing data constraints, but does not show the scatter in the growth histories among individual SMBHs. The blue dashed and dotted curves show the SMBH mass growth with a constant Eddington ratio of 1 and 5, respectively.

The high Eddington ratios at $z \gtrsim 6$ mean that SMBHs grow exponentially, although the e-folding time scale increases due to the decline in Eddington ratio with time. Consequently, the uncertainty in $M_{\bullet}$ increases as we move away from the starting redshift. This effect is the most significant for the 100 $M_{\odot}$ seed, whose growth was traced for more than 600 Myrs in one single direction. On the other hand, the average growth histories of GN-z11 and CEERS_1019 are only traced

<footnote>$^2$I.e., averaged over all active and dormant SMBHs of the same mass at the same time.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 10.** The *average* mass growth curves of GN-z11 (the blue solid curve and circle), CEERS_1019 (the red curve and circle), and a seed with $M_{\text{seed}} = 100 M_{\odot}$ at $z \sim 15$ (the orange curve and circle), predicted by TRINITY. The blue dashed and dotted curve shows the hypothetical mass growth histories of GN-z11 assuming constant Eddington ratios of 1 and 5, respectively. The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. These intervals reflect the ranges that the model predictions are allowed to vary, while keeping the underlying SMBH properties matching input data constraints within error bars (also see Section 2.2). The $z = 15$ progenitor masses of both GN-z11 and CEERS_1019 are $10^3-10^4 M_{\odot}$. See Section 5.2.6.</figcaption></figure>

for at most ~450 Myrs in each direction. At $z \sim 10.5$ and $z \sim 8.7$, the uncertainty in the $100 M_{\odot}$ seed's descendant mass marginally covers GN-z11 and CEERS_1019. This means that with current observational constraints on TRINITY, the potential progenitors of GN-z11, CEERS_1019 are unlikely to be as small as $100 M_{\odot}$ at $z \sim 15$.

*On average*, SMBHs grow at above (below) the Eddington rate above (below) $z \sim 14$. Thus, even though GN-z11 is accreting at five times the Eddington rate at the time of observation, such fast accretion episodes are likely intermittent (also see Section 5.2.1). This difference between the average and the instantaneous Eddington ratio is vital in estimating SMBH growth histories, given the exponential nature of early SMBH growth. As shown in Fig. 10, a constant Eddington ratio of 5 enables the fast mass buildup of ~$10^{10} M_{\odot}$ within ~120 Myrs. Conversely, it would take ~500 Myrs to achieve the same amount of mass growth with the Eddington rate. The drastically different SMBH growth histories also point to different SMBH progenitor masses, e.g. the constant super-Eddington growth implies the possibility of creating a GN-z11-like SMBH from with a $100 M_{\odot}$ progenitor at $z \sim 12$, which is strongly disfavored by the constant Eddington-limited scenario. When inferring potential SMBH progenitor masses, it is thus critical to use more accurate estimates of average Eddington ratios rather than instantaneous values. Since it is not possible to measure the average historical Eddington ratios for individual SMBH, models like TRINITY are essential in inferring average Eddington ratios for high-redshift SMBHs based on as many data constraints as possible.

We then discuss the potential seeding mechanisms for early SMBHs like GN-z11 and CEERS_1019. As can be seen from Fig. 10, the typical $z \sim 15$ progenitor masses of GN-z11 and CEERS_1019 are ~$10^3-10^4 M_{\odot}$, although CEERS_1019 could have grown from seeds as small as ~300 $M_{\odot}$ at $z \sim 15$. Such a progenitor mass range is consistent with the scenario where the seeds are created via runaway mergers in nuclear star clusters (see e.g. Devecchi et al. 2012; Lupi

et al. 2014; Kroupa et al. 2020; Schleicher et al. 2022; Kritos, Berti & Silk 2023, also see Topping et al. 2024 for observational evidence of such clusters at $z \sim 6$). In addition, when SMBH seeds form in star clusters, their dynamics are determined by the cluster mass rather than the seed mass, which makes them less susceptible to erratic motions than Population III star (PopIII) remnants. Compared to PopIII remnant seeds, they are thus also more likely to stay in high gas density regions and experience consistent growth as predicted by TRINITY. We thus believe that *if the SMBHs are seeded at z ~ 15*, the nuclear star cluster scenario is more plausible for objects like GN-z11 and CEERS_1019 than the PopIII remnant scenario. But since the 68th percentile range of GN-z11's and CEERS_1019's progenitor masses cover ~$10^4 M_{\odot}$ at $z \sim 11-14$, we *cannot* rule out the supermassive star scenario that would produce such heavy seed SMBHs (see e.g. Inayoshi, Visbal & Haiman 2020 and the references therein). While supermassive stars of ~$10^5 M_{\odot}$ formed via strong photodissociation radiation are predicted to have number densities below those of the detected high-z AGN, very massive/supermassive stars formed via dynamical heating are expected to be more numerous and have masses ~$10^4 M_{\odot}$ (see e.g. Wise et al. 2019; Regan et al. 2020; Lupi, Haiman & Volonteri 2021).

## 5.3 Constraints on the progenitor masses and growth rates of $M_{\bullet} \gtrsim 10^9 M_{\odot}$ SMBHs at $z \sim 9$

In this section, we compare the potential progenitor masses and growth rates of $M_{\bullet} \gtrsim 10^9 M_{\odot}$ SMBHs from Pacucci & Loeb (2022) and TRINITY. Given the difference between TRINITY and the Monte Carlo simulation by Pacucci & Loeb, we will only make qualitative statements on any differences to avoid overinterpretation. Pacucci & Loeb (2022) argued that the detection of $M_{\bullet} \gtrsim 10^9 M_{\odot}$ SMBHs at $z \gtrsim 9$ would strongly disfavour the PopIII remnant seeding mechanism, due to the insufficient time to grow substantially within the Eddington limit. Based on their Monte Carlo simulation, they conclude that the mean seed mass, $M_{\text{seed}}$, for $M_{\bullet} \gtrsim 10^9 M_{\odot}$ SMBHs at $z \gtrsim 9$ is $\log M_{\text{seed}} = 5.34$ at $z \sim 25$. Since the TRINITY model starts at $z \sim 15$, it is difficult to extrapolate towards higher redshifts with TRINITY. We thus opt to redo the Monte Carlo simulation as described in Pacucci & Loeb (2022), and compare their $z \sim 15$ progenitor masses of $M_{\bullet} \gtrsim 10^9 M_{\odot}$ SMBHs at $z \gtrsim 9$ with TRINITY predictions. By redoing this Monte Carlo simulation, we found that the $z \sim 15$ progenitor masses peak around $\log M_{\bullet, z \sim 15} \sim 7$. To calculate the $z \sim 15$ progenitor masses from TRINITY, we make use of the predicted mass-independence of the Eddington ratio at high redshifts. Specifically, the $z \sim 9$ descendant of GN-z11 is predicted to have $\log M_{\bullet} \sim 7$, according to Fig. 10. Its $z \sim 15$ progenitor has a mass of $\log M_{\bullet} \sim 4$. Given that Eddington ratio is nearly mass-independent, a $z \sim 15$ SMBH must be at least $\log M_{\bullet} \sim 4 + (9-7) = 6$ to reach $10^9 M_{\odot}$ by $z \sim 9$. This is roughly 1 dex lower than the typical $z \sim 15$ progenitor masses from Pacucci & Loeb (2022). Several factors could lead to this difference: (1) Pacucci & Loeb (2022) adopted a flat prior in $\log M_{\text{seed}}$, which is very different from the high-redshift BHMFs inferred by TRINITY (see Fig. 4). Consequently, the massive seeds will be given much larger weights in Pacucci & Loeb (2022) than in TRINITY; 2) In Pacucci & Loeb (2022), SMBHs are not allowed to accrete at super-Eddington rates at any time. On the contrary, SMBHs are predicted to accrete at super-Eddington rates at $z \gtrsim 14$ in TRINITY, even when averaged over active and dormant phases. This allows smaller SMBHs to grow into $\gtrsim 10^9 M_{\odot}$ SMBHs by $z \gtrsim 9$ in TRINITY.
With constant or slowly changing Eddington ratios, SMBH experience exponential growth:

$$M_{\bullet}(t) = M_{\text{seed}} \times \exp \left( \mathcal{D}\eta \frac{1-\epsilon}{\epsilon} \frac{\Delta t}{t_{\text{Edd}}} \right), \quad (11)$$

where $\mathcal{D}$ is the AGN duty cycle, i.e. the fraction of time during which SMBHs are accreting at an Eddington ratio $\eta$. $\epsilon$ is the AGN radiative efficiency, $\Delta t$ is the time interval between SMBH seeding and the time of observation $t$, and $t_{\text{Edd}} \equiv 450$ Myr. Since it is the product of $\mathcal{D}$, $\eta$, and $(1-\epsilon)/\epsilon$ that determines the growth rate, we define:

$$\lambda \equiv \mathcal{D}\eta \frac{1-\epsilon}{\epsilon}, \quad (12)$$

so that

$$M_{\bullet}(t) = M_{\text{seed}} \times \exp \left( \lambda \frac{\Delta t}{t_{\text{Edd}}} \right), \quad (13)$$

and compare the average $\lambda$ values from TRINITY with the peak value of the $\lambda$ distribution from the Monte Carlo simulation by Pacucci & Loeb (2022). Among $\gtrsim 10^9$ M$_{\odot}$ SMBHs by $z \sim 9$, the $\lambda$ distribution from Pacucci & Loeb (2022) peaks at $\lambda \sim 8.5$. In TRINITY, the average $\lambda$ decreases from $\sim 60$ at $z \sim 15$ to $\sim 8.8$ by $z \sim 9$. Hence, SMBHs in TRINITY grow much faster than simulated by Pacucci & Loeb (2022). Given the similar AGN efficiency and duty cycle values, this difference is mostly caused by disallowing super-Eddington accretion in Pacucci & Loeb (2022). This difference in growth rate directly affects the plausible progenitor masses for $\gtrsim 10^9$ M$_{\odot}$ SMBHs by $z \sim 9$ and could give different implications on SMBH seeding mechanisms. This comparison reemphasizes the importance of characterizing early SMBH growth rates in determining their seeding scenarios.

## 5.4 Constraining power of $z \gtrsim 6$ quasar luminosity functions

As mentioned in Section 2.1, there is no $z > 6.5$ SMBH data in the observational constraints for the current TRINITY model. Therefore, our predictions at $z \gtrsim 7$ are effectively extrapolations based on SMBH data at $z \le 6.5$ and galaxy data between $0 \le z \le 10$. Consequently, we see increasing statistical uncertainties towards higher redshifts in many predictions (e.g. Figs 3 and 4). However, new instruments like JWST, Euclid, and Roman Space Telescope will provide unprecedented amounts of $z > 6.5$ SMBH observations. Even before these new instruments, large ground-based telescopes already give constraints on QLFs up to $z \sim 7$ (e.g. Matsuoka et al. 2023). These high-redshift data will be vital for understanding SMBH formation and evolution in the early Universe. In this section, we demonstrate this by quantifying the constraining power of $z > 6.5$ quasar luminosity functions (QLFs) on the halo–galaxy–SMBH connection in TRINITY.

In the top panel of Fig. 11, we show the uncertainty in $M_{\bullet}$ as a function of $M_{\text{peak}}$ and $z$ from the fiducial TRINITY model, where no SMBH data at $z \gtrsim 6.5$ were included to constrain our model. The $M_{\bullet}$ uncertainty is a quadratic sum of: (1) the statistical uncertainty from the MCMC algorithm, which reflects the constraining power of all the data points given their values and uncertainties; and (2) the typical systematic uncertainty in $M_{\bullet}$ mass measurement from virial estimates, i.e. 0.5 dex. Due to the lack of high-redshift SMBH data, the $M_{\bullet}$ uncertainty increases significantly towards higher redshifts at $z \gtrsim 6$, reaching $\sim 1$ dex at $z \sim 10$ and $\gtrsim 2$ dex at $z \gtrsim 12$.

We then added mock $z = 6–10$ QLFs to quantify their constraining power on $z \gtrsim 6$ halo–galaxy–SMBH connection. These QLFs are extrapolated based on the modeling of $z \lesssim 6$ QLFs

<figure><img src="image_17.png" /><figcaption>**Figure 11.** Top Panel: The uncertainty in black hole mass ($M_{\bullet}$), as a function of $M_{\text{peak}}$ and $z$, with no quasar luminosity function (QLF) constraints between $z = 6–10$. Bottom Panel: same as the top panel, but with $z = 6–10$ mock QLF constraints from Roman High Latitude Wide Area Survey (~2000 deg$^2$). The uncertainties are the quadratic sums of the those from MCMC process and the typical $M_{\bullet}$ uncertainty from virial estimates, which is 0.5 dex. The MCMC uncertainties reflect the ranges that the model predictions are allowed to vary, while keeping the underlying SMBH properties matching input data constraints within error bars (also see Section 2.2). With the additional constraints from $z = 6–10$ mock QLF, the MCMC uncertainty in SMBH mass is greatly reduced between $z = 6–15$. See Section 5.4.</figcaption></figure>

by Shen et al. (2020). These extrapolated QLFs include both obscured and unobscured AGNs. We determine the upper and lower luminosity limits of the QLFs assuming the characteristics of Roman High Latitude Wide Area Survey (HLS, ~2000 deg$^2$); and (2) the Euclid Wide Survey (EWS, ~14500 deg$^2$). To calculate the minimum luminosities of the mock QLFs, we convolved the quasar SED from Temple, Hewett & Banerji (2021) with the filter transmission curves, and compare the fluxes with the near infrared magnitude limits of HLS (26.5 mag; Wang et al. 2022) and EWS (24.5 mag; Euclid Collaboration 2022). This yields minimum bolometric quasar luminosities of $10^{44.7}$ and $10^{45.5}$ erg s$^{-1}$ for HLS and EWS, respectively. We truncate the bright end of mock QLFs at the luminosities where only one quasar will be detected in the redshift-luminosity bins. This maximum luminosity value is further capped at $10^{48}$ erg s$^{-1}$, which is the Eddington luminosity for a $10^{10}$ M$_{\odot}$ SMBH. The maximum luminosity is $10^{48}$ erg s$^{-1}$ at $z = 6$, and $10^{46.5}$ and $10^{47.5}$ erg s$^{-1}$ at $z = 10$ for HLS and EWS, respectively.
As shown by the bottom panel of Fig. 11, the $z = 6–10$ QLFs from HLS will greatly reduce the $M_*$ uncertainty in the early Universe. In particular, the total $M_*$ uncertainties are now dominated by the virial estimate uncertainties in $M_*$ for haloes with $M_{peak} \gtrsim 10^{11}$ M$_\odot$ at $z \lesssim 10$. This means that $z = 6–10$ QLFs strongly constrain SMBH growth history between $z = 6–10$ when combined with all the other existing data. $z = 6–10$ SMBH growth histories are so strongly constrained that we can precisely calculate $z = 6–10$ SMBH masses by rewinding the growth histories of their $z \lesssim 6$ descendants. As such, our knowledge about $z = 6–10$ SMBH masses only depends on the uncertainties of their $z \lesssim 6$ descendant SMBH masses, i.e. the typical systematic uncertainty in $M_*$ mass measurement from virial estimates. We do not show the results with EWS, because there is only very slight improvement from HLS to EWS. This improvement comes from the better sky coverage of Euclid, which provides better statistics for bright quasars.

The inclusion of high-redshift QLFs does not significantly improve the constraints on early SMBH growth histories. In Fig. 12, we show the average SMBH mass growth curves for the same objects as in Fig. 10, but now with simulated $z = 6–10$ QLFs as additional constraints. The uncertainties in both progenitor and descendant masses of these objects only decrease slightly. This is because the redshift dependence of AGN Eddington ratio distributions is already well constrained (see Section 2.1), and the addition of QLFs ends up mostly constraining SMBH mass. But we caution that the high-redshift AGN Eddington ratios in TRINITY are effectively extrapolated from low-redshift data. The shaded regions shown in Figs 10 and 12 thus only represent the statistical uncertainties from the MCMC process, but do not include the systematic uncertainties from our input assumptions for such extrapolations, e.g. the exact parametrization of Eddington ratio distribution shape as a function of redshift.

# 6 CONCLUSIONS

In this work, we updated the TRINITY model of the halo–galaxy–SMBH connection, by adding the latest $9 \lesssim z \lesssim 13$ galaxy UV luminosity functions from JWST to our galaxy + SMBH data compilation. This additional constraint does not qualitatively change results presented in previous TRINITY papers (Zhang et al. 2023a, b, c), but does increase the number of galaxies and SMBHs at fixed masses beyond $z \sim 9$. With the updated TRINITY model, we predict the $z > 6 M_*$–$M_*$ relation for AGNs, the number density of massive ($M_* > 10^9 M_\odot$) SMBHs, and the host halo mass of massive SMBHs at $z \sim 6$. **Key results are as follows:**

(i) From $7 \lesssim z \lesssim 10$, the median $M_*$–$M_*$ relation for AGNs ($\log L_{bol} \ge 45$) is biased higher than the intrinsic relation for all SMBHs due to Lauer bias. The $M_*$–$M_*$ relation for AGNs is also roughly redshift-independent. Some of the $z \gtrsim 7$ AGNs detected by JWST (i.e. GN-z11 and CEERS-1019) are broadly consistent with our predicted $M_*$–$M_*$ relation given their AGN luminosities and mass measurement uncertainties, but UHZ1’s $M_*/M_*$ ratio is too high to match TRINITY predictions. (Section 4.1, Fig. 2).

(ii) The cumulative numbers of $M_* > 10^9$ and $10^{10} M_\odot$ black holes above a certain redshift increase by five orders of magnitude from $z \sim 10$ to $z \gtrsim 2$, due to fast SMBH growth. Massive SMBH numbers increase slowly below $z \sim 2$, because SMBH growth slows down and there is less comoving observable volume. Including observational scatter in measured $M_*$ inflates the numbers of massive black holes by a factor of $\sim 3–40$ due to Eddington bias. This boost is stronger

<figure><img src="image_18.png" /><figcaption>**Figure 12.** Similar average SMBH mass growth curves to Fig. 10, but with extra constraints from simulated $z = 6–10$ QLFs from Roman High Latitude Wide Area Survey. The shaded regions show the 68 per cent confidence intervals inferred from the model posterior distribution. These intervals reflect the ranges that the model predictions are allowed to vary, while keeping the underlying SMBH properties matching input data constraints within error bars (also see Section 2.2). The $z = 6–10$ mock QLFs only slightly improves our constraints on early SMBH growth histories. See Section 5.4.</figcaption></figure>

for more massive black holes at lower redshifts, where the SMBH mass functions are steeper. (Section 4.2, Figs 3–4).

(iii) $z \sim 6$ bright quasars with $M_* > 10^9$ and $10^{10} M_\odot$ are likely living in haloes with $M_{peak} \sim (2–3)$ and $(3–5) \times 10^{12} M_\odot$, which are $\sim 4.5–5.5\sigma$ peaks in Gaussian random fields. The $z=0$ relics of these quasars are typically hosted by haloes with $M_{peak} \gtrsim 3 \times 10^{14} M_\odot$. However, these host haloes only make up $\lesssim 1$ per cent of all the haloes at this mass scale, making it very difficult to find such relic SMBHs. We expect to find only $\lesssim 1$ such relic SMBH within $z \lesssim 0.05$ (Section 4.3, Fig. 5).

We also compared our predictions at $z \gtrsim 9$ with the two highest-redshift AGNs detected by JWST, GN-z11 ($z = 10.6$) and CEERS-1019 ($z = 8.7$). Our main findings are:

(i) The SMBH masses ($M_*$) and Eddington ratios ($\eta$) of GN-z11 and CEERS-1019 are consistent with TRINITY’s prediction that AGNs with $6 \lesssim \log M_* \lesssim 7$ and $\eta \gtrsim 1$ dominate the AGN luminosity function at $\log L_{bol} \sim 45$ and $z \sim 9–10$ (Section 5.2.1, Fig. 6).

(ii) Although both GN-z11 and CEERS-1019 are consistent with the predicted $M_*$–$M_*$ relation for luminosity-limited AGN samples, TRINITY predicts relatively small probabilities of detecting such objects in deep JWST surveys (0.02–0.04). In other words, there is a moderate tension between expected numbers of GN-z11 and CEERS-1019 analogues and the actual detections. However, the additional constraint from GN-z11 and CEERS-1019 does not change TRINITY modelling results significantly. Therefore, more observations are needed to better constrain TRINITY’s behaviour at highest redshifts.

(iii) The host halo masses are estimated to be $\sim (2–3) \times 10^{11} M_\odot$ and $\sim (3–5) \times 10^{11} M_\odot$ for GN-z11 and CEERS-1019, respectively. These halo masses correspond to $\gtrsim 7\sigma$ peaks in the density field at those redshifts, which are extremely biased regions of the Universe. The $z=0$ descendants of GN-z11 and CEERS-1019 are likely living in $(1–3) \times 10^{14} M_\odot$ haloes, but they are also very rare, similar to the local relics of $z=6$ quasars (Section 5.2.5, Fig. 9).

(iv) On average, SMBHs grow at above (below) the Eddington rate above (below) $z \sim 14$. The predicted average growth histories
make GN-z11 a potential progenitor of $z \sim 6$ quasars with $\log M_\odot \gtrsim 9$, which is less likely to be the case for CEERS_1019.

(v) The fraction of galaxies hosting AGNs with bolometric luminosities above $10^{45}$ erg s$^{-1}$, $f_{\text{gal}}$ is a strong function of galaxy mass. At $z \gtrsim 6$, TRINITY predicts that $f_{\text{gal}}$ increases from $\lesssim 10^{-6}$ at $M_\star \sim 10^9 M_\odot$ to $\sim 1$ at $M_\star \sim 3 \times 10^{10} M_\odot$ (Section 5.2.3, Fig. 8).

(vi) The cumulative $z=8-12$ AGN luminosity function calculated from CEERS_1019 and GN-z11 does not significantly change TRINITY predictions. To impose a strong enough constraint on TRINITY, more detections of AGNs above $L_{\text{bol}} = 10^{45}$ erg s$^{-1}$ and/or better measurements of host galaxy mass would be needed. Alternatively, non-detection of such AGNs over a larger area would be necessary to confirm current TRINITY predictions.

The comparison between TRINITY predictions and GN-z11 and CEERS_1019 implies that runaway mergers in nuclear star clusters or very massive stars mediated by dynamical heating would be a viable SMBH seeding mechanism, assuming that these seeds formed at $z \sim 15$. (Section 5.2.6, Fig. 10)

Finally, we showed that future $6 < z < 10$ AGN data will provide significant constraints on the early halo-galaxy-SMBH mass connection. Specifically, $z=6-10$ quasar luminosity functions (QLFs) measured by future wide area surveys (e.g. the Roman High Latitude Wide Area Survey and the Euclid Wide Survey) will reduce SMBH mass uncertainty from $\sim 1.0$ to $\sim 0.5$ dex for haloes with $M_{\text{peak}} \gtrsim 10^{11} M_\odot$ at $z \lesssim 10$. (Section 5.4, Fig. 11).

## ACKNOWLEDGEMENTS

We thank the anonymous referee for their constructive comments on the manuscript. We thank Gurtina Besla, Frank van den Bosch, Haley Bowden, Alison Coil, Anna-Christina Eilers, Ryan Endsley, Sandy Faber, Nickolay Gnedin, Richard Green, Jenny Greene, Kate Grier, Melanie Habouzit, Kevin Hainline, Andrew Hearin, Yun-Hsin Huang, Raphael Hviding, Stephanie Juneau, David Koo, Andrey Kravtsov, Daniel Lawther, Jianwei Lyu, Chung-Pei Ma, Daniel Mortlock, Vasileios Paschalidis, Joel Primack, Eliot Quataert, George Rieke, Marcia Rieke, Jan-Torge Schindler, Xuejian Shen, Yue Shen, Dongdong Shi, Rachel Somerville, Fengwu Sun, Yoshihiro Ueda, Marianne Vestergaard, Ben Weiner, Christina Williams, Charity Woodrum, and Minghao Yue for very valuable discussions.

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).

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de).

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 AVAILABILITY

The parallel implementation of TRINITY, the compiled data sets, and the posterior distribution of model parameters are available at https://github.com/HaowenZhang/TRINITY.

## REFERENCES

Agarwal B., Davis A. J., Khochfar S., Natarajan P., Dunlop J. S., 2013, MNRAS, 432, 3438

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 A Rev., 56, 93

Angulo R. E., Springel V., White S. D. M., Cole S., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 425, 2722

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

Bellovary J. M., Cleary C. E., Munshi F., Tremmel M., Christensen C. R., Brooks A., Quinn T. R., 2019, MNRAS, 482, 2913

Blandford R. D., McKee C. F., 1982, ApJ, 255, 419

Bogdan A. et al., 2023, Nat. Astron., 8, 126

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&A Rev., 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., Storchi-Bergmann T., 2000, ApJ, 533, 682

Chabrier G., 2003, PASP, 115, 763

Chiang Y.-K., Overzier R., Gebhardt K., 2013, ApJ, 779, 127

Choksi N., Behroozi P., Volonteri M., Schneider R., Ma C.-P., Silk J., Moster B., 2017, MNRAS, 472, 1526

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

Dekel A., Lapiner S., Dubois Y., 2019, preprint (arXiv:1904.08431)

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

Devecchi B., Volonteri M., Rossi E. M., Colpi M., Portegies Zwart S., 2012, MNRAS, 421, 1465

Diemer B., 2018, ApJS, 239, 35

Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662

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

Eddington A. S., 1913, MNRAS, 73, 359

Eisenstein D. J., Hu W., 1998, ApJ, 496, 605

Euclid Collaboration, 2019, A&A, 631, A85

Euclid Collaboration, 2022, A&A, 662, A112

Fan X., 1999, AJ, 117, 2528

Fanidakis N., Macciò A. V., Baugh C. M., Lacey C. G., Frenk C. S., 2013, MNRAS, 436, 315

Feng Y., Di Matteo T., Croft R., Khandai N., 2014, MNRAS, 440, 1865

Ferrarese L., Ford H., 2005, Space Sci. Rev., 116, 523

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

Finkelstein S. L. et al., 2019, ApJ, 879, 36

<footnote n="3">³https://apps.automeris.io/wpd/</footnote>
Gebhardt K. et al., 2000, ApJ, 539, L13

Goulding A. D. et al., 2023, ApJ, 955, L24

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

Harikane Y. et al., 2023a, ApJ, 959, 39

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

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

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

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

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

Inayoshi K., Visbal E., Haiman Z., 2020, ARA&A, 58, 27

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

Izumi T. et al., 2021, ApJ, 914, 36

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

Jiang L. et al., 2016, ApJ, 833, 222

Kashino D., Lilly S. J., Matthee J., Eilers A.-C., Mackenzie R., Bordoloi R., Simcoe R. A., 2023, ApJ, 950, 66

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

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

Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340

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

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

Kritos K., Berti E., Silk J., 2023, Phys. Rev. D, 108, 083012

Kroupa P., Subr L., Jerabkova T., Wang L., 2020, MNRAS, 498, 5652

Kulkarni G., Worseck G., Hennawi J. F., 2019, MNRAS, 488, 1035

Labbé I. et al., 2023, Nature, 616, 266

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

Larson R. L. et al., 2023, ApJ, 953, L29

Latif M. A., Whalen D. J., Khochfar S., Herrington N. P., Woods T. E., 2022, Nature, 607, 48

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

Li Y. et al., 2007, ApJ, 665, 187

Li J. et al., 2024, preprint (arXiv:2403.00074)

Lupi A., Colpi M., Devecchi B., Galanti G., Volonteri M., 2014, MNRAS, 442, 3616

Lupi A., Haiman Z., Volonteri M., 2021, MNRAS, 503, 5046

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

Maiolino R. et al., 2023a, Nature, 627, 59

Maiolino R. et al., 2023b, preprint (arXiv:2308.01230)

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

Matsuoka Y. et al., 2023, ApJ, 949, L42

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

McGree I. D. et al., 2013, ApJ, 768, 105

McGree I. D., Fan X., Jiang L., Cai Z., 2018, AJ, 155, 131

McLure R. J., Dunlop J. S., 2004, MNRAS, 352, 1390

McLure R. J., Jarvis M. J., 2002, MNRAS, 337, 109

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

Merloni A., 2004, MNRAS, 353, 1035

Merloni A. et al., 2014, MNRAS, 437, 3550

Natarajan P., Pacucci F., Ferrara A., Agarwal B., Ricarte A., Zackrisson E., Cappelluti N., 2017, ApJ, 838, 117

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

Pacucci F., Loeb A., 2022, MNRAS, 509, 1885

Peterson B. M., 1993, PASP, 105, 247

Pillepich A. et al., 2018, MNRAS, 475, 648

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

Regan J. A., Wise J. H., Woods T. E., Downes T. P., O’Shea B. W., Norman M. L., 2020, Open J. Astrophys., 3, 15

Richards G. T. et al., 2002, AJ, 123, 2945

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

Ross N. P. et al., 2012, ApJS, 199, 3

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

Schleicher D. R. G. et al., 2022, MNRAS, 512, 6192

Schneider R., Valiante R., Trinca A., Graziani L., Volonteri M., 2023, MNRAS, 526, 3250

Scholtz J. et al., 2023, preprint (arXiv:2306.09142)

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. et al., 2019, ApJ, 873, 35

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

Shimasaku K., Izumi T., 2019, ApJ, 872, L29

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, Ap&A, 331, L1

Silk J., Begelman M. C., Norman C., Nusser A., Wyse R. F. G., 2024, ApJ, 961, L39

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

Soltan A., 1982, MNRAS, 200, 115

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

Springel V. et al., 2005, Nature, 435, 629

Tacchella S. et al., 2023, ApJ, 952, 74

Tee W. L., Fan X., Wang F., Yang J., Malhotra S., Rhoads J. E., 2023, ApJ, 956, 52

Temple M. J., Hewett P. C., Banerji M., 2021, MNRAS, 508, 737

Tenneti A., Di Matteo T., Croft R., Garcia T., Feng Y., 2018, MNRAS, 474, 597

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

Topping M. W. et al., 2024, MNRAS, 529, 3301

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

Tremmel M., Karcher M., Governato F., Volonteri M., Quinn T. R., Pontzen A., Anderson L., Bellovary J., 2017, MNRAS, 470, 1121

Trinca A., Schneider R., Valiante R., Graziani L., Zappacosta L., Shankar F., 2022, MNRAS, 511, 616

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

Vestergaard M., Osmer P. S., 2009, ApJ, 699, 800

Vestergaard M., Peterson B. M., 2006, ApJ, 641, 689

Volonteri M., Stark D. P., 2011, MNRAS, 417, 2085

Volonteri M., Reines A. E., Atek H., Stark D. P., Trebitsch M., 2017, ApJ, 849, 155

Volonteri M., Habouzit M., Colpi M., 2023, MNRAS, 521, 241

Wang F. et al., 2019, ApJ, 884, 30

Wang F. et al., 2021, ApJ, 907, L1

Wang Y. et al., 2022, ApJ, 928, 1

Wang F. et al., 2023, ApJ, 951, L4

Wise J. H., Regan J. A., O’Shea B. W., Norman M. L., Downes T. P., Xu H., 2019, Nature, 566, 85

Yang J. et al., 2016, ApJ, 829, 33

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

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

Zhang H., Behroozi P., Volonteri M., Silk J., Fan X., Aird J., Yang J., Hopkins P. F., 2023b, MNRAS, 529, 2777

Zhang H., Behroozi P., Volonteri M., Silk J., Fan X., Hopkins P. F., Yang J., Aird J., 2023c, MNRAS, 518, 2123

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