Trinity IV: Predictions for Supermassive Black Holes at $z\gtrsim 7$

As an empirical model, Trinity infers the range of halo–galaxy–SMBH connections that match our unique compilation of galaxy$+$AGN data. Many predictions in this work come from the following two aspects of this connection:

Firstly, the intrinsic $M_{\bullet}$–$M_{*}$ relation evolves only mildly with redshift (see §4.1). The largest change occurs for low-mass galaxies (i.e., $M_{*}\sim 5\times 10^{9}M_{\odot}$) which have lower SMBH masses at higher redshifts; almost no change in SMBH masses with redshift occurs for larger 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_{\mathrm{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 that can simultaneously reproduce: 1) the probability distribution of AGN luminosity as a function of galaxy mass and redshift, i.e., QPDFs; 2) QLFs; and 3) the local $M_{\bullet}$–$M_{\mathrm{bulge}}$ relation. For example, a steeper $M_{\bullet}$–$M_{*}$ relation will not only violate the constraint from the local $M_{\bullet}$–$M_{\mathrm{bulge}}$ relation, but also induce faster growth of SMBHs at fixed galaxy mass, and thus overproduce AGN activity compared to observed QPDFs.

Secondly, the AGN Eddington ratio distribution becomes narrower towards higher redshifts, along 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_{\mathrm{bulge}}$ relation, lower Eddington ratios and/or broader Eddington ratio distributions will underpredict the number of SMBHs above certain luminosity limits, i.e., ABHMFs, at higher redshifts.

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 relatively non-evolving SMBH–galaxy relationship at those redshifts. Hence, a failure of Trinity to predict $z>6$ observations would indicate different evolution in the SMBH–galaxy mass relationship than that seen at $z<6$. As we will show in §4.1 and §5.2, the SMBH properties of CEERS1019 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 very 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.

Trinity makes self-consisitent and statistical halo–galaxy–SMBH connections from $z=0-10$. Instead of following individual haloes across cosmic time like the Universe Machine (Behroozi et al., 2019), Trinity is built upon halo statistics from N-body simulations. The halo–galaxy connection is made via the $\mathrm{SFR}$–$M_{\mathrm{peak}}$ relation, where $M_{\mathrm{peak}}$ is the peak halo mass of haloes throughout their assembly histories. With a given redshift-dependent $\mathrm{SFR}$–$M_{\mathrm{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_{\mathrm{bulge}}$) with a redshift-dependent but fixed $M_{\mathrm{bulge}}$–$M_{*}$ based on SDSS (Mendel et al., 2014) and CANDELS (Lang et al., 2014) observations. We parametrize the $M_{\bullet}$–$M_{\mathrm{bulge}}$ relation as a redshift-dependent power-law for all galaxies, i.e., SMBH hosts and non-hosts. The log-normal scatter around the $M_{\bullet}$–$M_{\mathrm{bulge}}$ relation, $\sigma_{\mathrm{BH}}$, is a free parameter to be constrained in MCMC. We assume that $\sigma_{\mathrm{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 parameterize the $M_{\bullet}$–$M_{*}$ relation as a power-law. SMBH occupation fractions are parameterized as a sigmoid function of halo mass, which evolves with redshift. Together, these choices fully parameterize the SMBH mass distribution as a function of galaxy mass and redshift ($P(M_{\bullet}|M_{\ast},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_{\mathrm{bol}}|M_{\bullet},M_{\ast},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 et al. 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. The whole parametrization of Trinity contains 54 free parameters.

In this work, we present various statistics of SMBHs at $z\gtrsim 6$ 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_{\mathrm{bol,lim}}$), we firstly calculate the SMBH mass function at each given stellar mass and redshift:

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

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

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

where $P\left(M_{\bullet}|M_{\rm peak},z\right)$ 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_{\mathrm{peak}}$, $P(M_{*}|M_{\mathrm{peak}},z)$:

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

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

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

As noted in §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_{\mathrm{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_{\mathrm{peak}}$ instead of the present day halo mass. These adjustments were based on the Bolshoi & Consuelo simulations (Klypin et al., 2011). We refer readers to Appendix G of Behroozi et al. (2013) for full details. With these calibrations, the halo statistics used in this work are suitable for studying the evolution of halos from $10^{10}M_{\odot}$ 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 assumes Planck cosmology (Planck Collaboration et al., 2016), and has a simulation box size of 1 Gpc$/h$. The mass of each dark matter particle is $1.51\times 10^{9}M_{\odot}/h$, which enables the simulation to resolve haloes down to $M_{\mathrm{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 growth histories of haloes hosting SMBHs with $M_{\bullet}>10^{9}M_{\odot}$ at $z\sim 6$, as well as the descendants of typical haloes hosting GNz-11 and CEERS_1019.

To constrain the 54 free parameters in Trinity, we have originally compiled a diverse dataset 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\leq z\leq 10$. The SMBH observables are: quasar luminosity functions (QLFs), quasar probability distribution functions (QPDFs), active black hole mass functions (ABHMFs), the $z=0$ $M_{\bullet}$–$M_{\mathrm{bulge}}$ relation, and the observed $M_{\bullet}$ distribution of high redshift bright quasars. These SMBH observables collectively cover a redshift range of $0\leq z\leq 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 previous 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).

The addition of $9\lesssim z\lesssim 13$ galaxy UVLFs does not induce any qualitative changes in the results presented by Zhang et al. (2023c), Zhang et al. (2023a), or Zhang et al. (2023b). Quantitatively, it does increase the number of UV-bright galaxies at $z\gtrsim 9$, and ultimately, the galaxy masses at fixed halo masses. This increase in galaxy number density and typical galaxy mass also leads to more SMBHs at fixed SMBH masses, because the redshift-dependent $M_{\bullet}$–$M_{*}$ relation is unchanged with the addition of new galaxy UVLFs. As shown in Fig. 1, this change in galaxy number density increases towards higher redshifts, especially above $z=10$. This is because the new JWST UVLFs (solid circles) are largely consistent with existing UVLFs in our data compilation (solid stars), but we lacked constraints at $z\sim 12$ in the original Trinity model. Without $z\sim 12$ constraints, the original Trinity underproduced UV-bright galaxies, leading to the larger difference at this redshift. 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 here.

In Zhang et al. (2023a), we showed the median $M_{\bullet}$–$M_{*}$ relation for bright quasars at $z\sim 6$, and concluded that overmassive SMBHs compared to the observed local $M_{\bullet}$–$M_{*}$ relation are consistent with arising purely from selection bias against less massive and thus fainter quasars. In Fig. 2, we extend our predictions to even higher redshifts, i.e., $6\lesssim z\lesssim 11$. The solid lines denote the intrinsic (i.e., including active and inactive SMBHs) median $M_{\bullet}$–$M_{*}$ relation at different redshifts. The non-solid curves are the $M_{\bullet}$–$M_{*}$ relations for quasars above various bolometric luminosity limits, which are labeled next to each line. Given their weak redshift dependence, we only show these quasar $M_{\bullet}$–$M_{*}$ relations at $z=11$ for visual clarity. The shaded regions correspond to $1-\sigma$ log-normal scatter of $\sim 0.6$ dex around these quasar $M_{\bullet}$–$M_{*}$ relations, which includes: 1) the intrinsic scatter around the $M_{\bullet}$–$M_{*}$ relation ($\sim 0.3$ dex, according to the best-fit 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_{\mathrm{bol}}>48$ and the deepest samples with $\log L_{\mathrm{bol}}>45$, also for clarity.

Similar to those at $z\sim 6$, the $M_{\bullet}$–$M_{*}$ 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 §2.1). When selected by bolometric luminosity, more massive SMBHs are brighter and naturally over-represented in the sample. At fixed bolometric luminosity limits, the $M_{\bullet}$–$M_{*}$ relations for quasars barely evolve with redshift, despite the evolution in the intrinsic $M_{\bullet}$–$M_{*}$ relation for all SMBHs. This is because at fixed SMBH mass, the Eddington ratio distribution evolves very mildly beyond $z\sim 6$. Thus, selecting quasars by their luminosities yields similar SMBH masses across different redshifts, which are more and more overmassive SMBHs compared to the lower intrinsic $M_{\bullet}$–$M_{*}$ relations towards higher redshifts.

We also note that the typical $M_{\bullet}$ for quasars depends only weakly on host galaxy mass below $\log M_{*}\sim 9.5$. This is because, on average (i.e., including their active and dormant phases), high-redshift SMBHs grow at slightly below the Eddington rate, and hyper-Eddington AGNs are very rare. This upper limit in Eddington ratios translates into a lower limit in the typical $M_{\bullet}$, regardless of host galaxy mass. For brighter quasars (e.g., $\log L_{\mathrm{bol}}\geq 48$), the typical SMBH mass remains above $3\times 10^{8}M_{\odot}$ even for galaxies below the same mass. This implies that if we detected any $\log L_{\mathrm{bol}}\geq 48$ quasars in a $M_{*}<3\times 10^{8}M_{\odot}$ galaxy, they either: 1) have Eddington ratios much higher than predicted by Trinity; and/or 2) they are in the outsized black hole galaxy (OBG) phase (Agarwal et al., 2013; Natarajan et al., 2017).

Of note, when quantifying this luminosity dependence of the $M_{\bullet}$–$M_{*}$ relation (i.e., the Lauer bias), we assume that every AGN above those luminosity limits is included in the samples, with SMBH mass measurements via virial estimates. In other words, we do not account for the possibility that these AGNs may be hidden in the light from host galaxies. This assumption is more likely to hold for brighter quasars. For fainter AGNs, e.g., $\log L_{\mathrm{bol}}\sim 45$, stronger contamination from the host galaxies makes it more difficult to measure SMBH masses, even if the AGN luminosities are technically above the limits of current instruments. As a result, overmassive SMBHs would still be preferred among fainter AGNs (see, e.g., Volonteri et al. 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. Therefore, our predicted $M_{\bullet}$–$M_{*}$ relations for fainter AGNs (e.g., with $\log L_{\mathrm{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), CEERS1019 (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 X-ray spectral information, modeling jointly luminosity and obscuring column density, and further assuming Eddington-limited accretion. Given that SMBHs accreting at five times the Eddington rate have been detected at similar redshifts (i.e., GNz-11), we adopt a higher $M_{\bullet}$ uncertainty of 0.7 dex (i.e., a factor of five) to reflect the uncertainty in Eddington ratio. Galaxy mass uncertainties are directly taken from these source papers. We note that the stellar masses have in several cases been determined without including an AGN component in the fit, and that the uncertainties in the star formation histories can give 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_{\mathrm{bol}}\sim 47$ compiled by Izumi et al. (2021), all these JWST AGNs have lower SMBH masses. This is consistent with Trinity’s prediction for the luminosity-dependent bias (i.e., the Lauer bias) of 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 compare our predictions with CEERS_00717 after more precise measurements of galaxy mass are available. We do note that both CEERS1019 and GN-z11 lie above our $M_{\bullet}$–$M_{*}$ relation for AGNs with $\log L_{\mathrm{bol}}>45$, at $\sim 1.5\sigma$ and $\sim 0.8\sigma$, respectively. This could be due to: 1) small number statistics; 2) the aforementioned selection bias due to smaller SMBHs being harder to detect against the host galaxy’s light; 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 §5.2, we compare CEERS1019 and GN-z11 with other $6\lesssim z\lesssim 11$ predictions from Trinity.

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_{\mathrm{bol}}\geq 45.7$, the dash-dotted line). But, at $z\gtrsim 10$, another super-Eddington AGN has been found, i.e., GN-z11. If UHZ1 is in fact 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 thus $M_{\bullet}$ as well) is only a $1\sigma$ lower limit from X-ray spectral fitting. Such values were reported instead of the actual best-fitting values due to the strong degeneracy between the intrinsic X-ray luminosity and gas column density used in the fitting process. Thus, choosing higher intrinsic luminosity values will push UHZ1 further away from Trinity’s prediction. This would imply that UHZ1 may be experiencing an OBG phase, where the SMBH is not evolving in sync with the host galaxy and could become more massive than the galaxy. If confirmed, such AGNs in OBGs will provide critical information for future versions of Trinity so that it can include these extreme objects.

The top panel of Fig. 3 shows the redshift density of very massive black holes ($M_{\bullet}>10^{9}M_{\odot}$ and $10^{10}M_{\odot}$), i.e., $dN/dz$, 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 black hole accretion rates. Below $z\sim 1-2$, $dN/dz$ decreases due to the combination of: 1) the decreasing observable volume; 2) the decline in SMBH growth. The middle panel of Fig. 3 shows cumulative numbers of very massive black holes at different redshifts in the observable Universe. 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), because there are more low-mass black holes that can be upscattered than high-mass black holes that can be downscattered. This elevates the observed number of massive black holes above the true number, as shown by the dashed curves in the top and middle panels of Fig. 3. Overall, an observational scatter of $0.5$ dex elevates observed massive black hole number counts 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, as is shown in 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},10^{10}M_{\odot}$ SMBHs increase by six 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$ (the bottom right panel of Fig. 2). The rarity of such SMBHs is also shown in the bottom panel of Fig. 3: we expect to see only $\lesssim 10$ 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). The detection of massive BHs at $z\gtrsim 10$ will therefore provide stringent constraints on SMBH seeding mechanisms, as well as 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 §5.3.

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 bolometric luminosity limits as a function of redshift, we adopted the magnitude limits from Ivezić et al. (2019) and Euclid Collaboration et al. (2022). We further convert observed magnitude limits in different wavebands (labeled in the bottom panel of Fig. 3) into bolometric luminosities, by assuming: 1) that the bolometric correction to the rest-frame $i$-band luminosity is 12 (Richards et al., 2006); and 2) the UV-to-optical quasar SED is a power-law with a power-law index of $\alpha_{\nu}=-0.44$ (Vanden Berk et al., 2001). The luminosity limit as a function of redshift is shown in the bottom panel of Fig. 3, and the differential (cumulative) number of 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,\mathrm{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 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 these high luminosity limits ($\gtrsim 10^{46}$ erg/s), we expect $\sim 40\%$($80\%$) of such SMBHs to be unobscured, based on the obscured fractions as functions of AGN X-ray luminosity given by Merloni et al. (2014) and Ueda et al. (2014). Below $z\sim 4$, massive black holes experience a significant decline in activity level. Consequently, smaller and smaller fractions of massive black holes at lower redshifts will lie above the flux limit of LSST.

We show the distribution of host halo masses for SMBHs with $M_{\bullet}>10^{9}M_{\odot}$ (solid curves) and $M_{\bullet}>10^{10}M_{\odot}$ (dashed curves) at $z=6$, in the top panel of Fig. 5. When using virial estimates to measure $M_{\bullet}$, there is random scatter in the measured $M_{\bullet}$ around the intrinsic values. So, we also show the host halo mass functions for observed $M_{\bullet}>10^{9}M_{\odot}$ in blue and red curves, respectively. In this comparison, we assume a random scatter of $0.5$ dex in measured $M_{\bullet}$. This random scatter induces Eddington bias, i.e., more lower-mass haloes hosting intrinsically smaller SMBHs have their observed $M_{\bullet}$ upscattered than the opposite. Thus, the host halo mass function (HMF) is both enhanced in normalization and shifted towards lower masses. The median host halo mass for SMBHs with intrinsic $M_{\bullet}>10^{9}M_{\odot}$ is $\sim 3\times 10^{12}M_{\odot}$, compared to $\sim 2\times 10^{12}M_{\odot}$ for SMBHs with observed $M_{\bullet}>10^{9}M_{\odot}$. For $M_{\bullet}>10^{10}M_{\odot}$ SMBHs, their host haloes are typically $\sim 5\times 10^{12}M_{\odot}$ and $\sim 3\times 10^{12}M_{\odot}$, for intrinsic and observed $M_{\bullet}$, respectively. But overall, host haloes of $M_{\bullet}>10^{9}(10^{10})M_{\odot}$ SMBHs only make up $\sim 10\%(1\%)$ of all the haloes (orange solid curve, calculated from the MDPL2 simulation) at their respective median host halo masses, which correspond to $\sim 4.5-5.5\sigma$ peaks in Gaussian random fields (Eisenstein & Hu, 1998; Diemer, 2018). In such haloes, gas reservoirs available for rapid SMBH growth could include 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 catalogs 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 studies investigating where their descendant SMBHs can be found in the local Universe (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 halos with the progenitor’s mass at $z=6$ that would host a $>10^{9}/10^{10}M_{\odot}$ SMBH. As shown in the bottom panel of Fig. 5, the descendant HMFs are much broader than the progenitor HMFs. This is due to the diversity in the assembly histories of individual haloes. At $z=0$, descendants with $M_{\bullet}(z=6)>10^{9}M_{\odot}$ typically have $M_{\mathrm{peak}}\sim 3-7\times 10^{14}M_{\odot}$. However, these haloes only make up $\lesssim 1\%$ of the haloes at similar mass scales (the orange solid line). Combined with the small comoving volume in the local Universe, we would expect only $\sim 1$ such relic SMBH closer than $z=0.05$. For $M_{\bullet}(z=6)>10^{10}M_{\odot}$, their host halo masses are similar to those with $M_{\bullet}(z=6)>10^{9}M_{\odot}$, but they are much rarer (i.e., a factor of $\sim 30-100$). This makes it even harder to find the descendant haloes of these more massive high-redshift SMBHs in the local Universe.

According to Trinity, the intrinsic $M_{\bullet}$–$M_{*}$ relation between $7\lesssim z\lesssim 10$ is relatively static, with the largest difference being that the relation is mildy steeper at higher redshifts than at lower redshifts (see also Volonteri & Stark 2011). This means that SMBHs are less massive than in local galaxies at fixed stellar masses below $5\times 10^{10}M_{\odot}$ (§4.1). Physically, this redshift evolution implies that SMBH growth in $M_{*}<10^{11}M_{\odot}$ galaxies is modestly slower than the host galaxy’s growth. 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., IllustrisTNG, 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 et al. 2017, Bellovary et al. 2019, and Dekel et al. 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\gtrsim 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).

In this section, we compare the potential progenitor masses and growth rates of $\gtrsim 10^{9}M_{\odot}$ SMBHs from Pacucci & Loeb (2022) and Trinity. Since the methodologies of Trinity and the Monte Carlo simulation by Pacucci & Loeb are very different, we will only make qualitative statements on any differences, to avoid over-interpretation. Pacucci & Loeb (2022) argued that the detection of $\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_{\mathrm{seed}}$, for $\geq 10^{9}M_{\odot}$ SMBHs at $z\gtrsim 9$ is $\log M_{\mathrm{seed}}=5.34$ at $z\sim 25$. Since the Trinity model starts at $z\sim 15$, it is impossible to extrapolate towards higher redshifts with Trinity. We thus opt to redo the Monte Carlo simulation as described in Pacucci & Loeb (2022), and compare the $z\sim 15$ progenitor masses of $\gtrsim 10^{9}M_{\odot}$ SMBHs at $z\gtrsim 9$ with Trinity.

We find after this exercise that the $z\sim 15$ progenitor masses peak around $\log M_{\mathrm{\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.5$, according to Fig. 9. Its $z\sim 15$ progenitor has a mass of $\log M_{\bullet}\sim 3$. Given that Eddington ratio is nearly mass-independent, a $z\sim 15$ SMBH must be at least $\log M_{\bullet}\sim 3+(9-7.5)=4.5$ to reach $10^{9}M_{\odot}$ by $z\sim 9$. This is roughly 2.5 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_{\mathrm{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 10.5$ 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:

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_{\mathrm{Edd}}\equiv 450$ Myr. Since it is the product of $\mathcal{D}$, $\eta$, and $(1-\epsilon)/\epsilon$ that determines the growth rate, we define:

so that

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$. In other words, SMBHs in Trinity grow significantly faster than simulated by Pacucci & Loeb (2022). Given the similar AGN efficiency and duty cycle values, this difference is mostly caused by allowing/disallowing super-Eddington accretion in Trinity/Pacucci & Loeb (2022). As mentioned in the last paragraph, 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. Our comparison here demonstrates again the importance of characterizing early SMBH growth rates in determining their seeding scenarios.

As mentioned in §2.1, there is currently no $z>6.5$ SMBH data in the observational constraints for Trinity. Therefore, our predictions at $z\gtrsim 7$ are effectively extrapolations based on SMBH data at $z\leq 6.5$, in addition to galaxy data between $0<z<10$. Consequently, we see increasing statistical uncertainties towards higher redshifts in many predictions (e.g., Figs. 3 and 4). However, new instruments like JWST, Roman Space Telescope, and Euclid will provide an unprecedented amount 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_{\mathrm{peak}}$ and $z$ from the fiducial Trinity model, where no SMBH data at $z\gtrsim 6$ 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 when their values and uncertainties are perfectly accurate; 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 see whether/how much these data can help constrain the halo–galaxy–SMBH connection at $z\gtrsim 6$. These QLFs are extrapolated based on the modeling (global fit B) of $z\lesssim 6$ QLFs by Shen et al. (2020). These extrapolated QLFs are based on obscuration-corrected lower-redshift QLFs, so in principle 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, $\sim 2000$ deg²); and 2) the Euclid Wide Survey (EWS, $\sim 14500$ deg²). To calculate the minimum luminosities of the mock QLFs, we convolved the quasar SED from Temple et al. (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 et al. 2022). This yields minimum bolometric quasar luminosities of $10^{44.7}$ and $10^{45.5}$ erg/s 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, which is the Eddington luminosity for a $10^{10}M_{\odot}$ SMBH. For HLS (EWS), the maximum luminosity is $10^{48}$ ($10^{48}$) erg/s at $z=6$, and $10^{46.5}$ ($10^{47.5}$) erg/s at $z=10$.

As shown by the bottom panel of Fig. 11, the $z=6-10$ QLFs from HLS will greatly reduce the $M_{\bullet}$ uncertainty in the early Universe. In particular, the total $M_{\bullet}$ uncertainties are now dominated by the virial estimate uncertainties in $M_{\bullet}$ for haloes with $M_{\mathrm{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 they are combined with (the high-redshift extrapolation of) all the other existing data included by Trinity. $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 accuracy of their $z\lesssim 6$ descendant SMBH masses. 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 also 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. 9, 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 §2.1), and the addition of QLFs ends up mostly constraining SMBH mass. Nonetheless, we do caution that the high-redshift AGN Eddington ratios in Trinity are effectively extrapolated from low-redshift data. The shaded regions shown in Figs. 9 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 parameterization of Eddington ratio distribution shape as a function of redshift.