# Tip of the Iceberg: Overmassive Black Holes at $4 < z < 7$ Found by JWST Are Not Inconsistent with the Local Relation

Junyao Li, John D. Silverman, Yue Shen, Marta Volonteri, Knud Jahnke, Ming-Yang Zhuang, Matthew T. Scoggins, Xuheng Ding, Yuichi Harikane, Masafusa Onoue, et al.

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

## To cite this version:

Junyao Li, John D. Silverman, Yue Shen, Marta Volonteri, Knud Jahnke, et al.

HAL Id: insu-05101687

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

Submitted on 7 Jun 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

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

# Tip of the Iceberg: Overmassive Black Holes at $4 < z < 7$ Found by JWST Are Not Inconsistent with the Local $M_{BH}-M_{\star}$ Relation

Junyao Li¹, John D. Silverman²,³,⁴,⁵, Yue Shen¹,⁶, Marta Volonteri⁷, Knud Jahnke⁸, Ming-Yang Zhuang¹, Matthew T. Scoggins⁹, Xuheng Ding¹⁰, Yuichi Harikane¹¹, Masafusa Onoue²,⁴,¹², and Takumi S. Tanaka²,³,⁴

¹Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA; junyaoli@illinois.edu

²Kavli Institute for the Physics and Mathematics of the Universe, The University of Tokyo, Kashiwa, 277-8583 (Kavli IPMU, WPI), Japan

³Department of Astronomy, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan

⁴Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

⁵Center for Astrophysical Sciences, Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA

⁶National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA

⁷Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis bd Arago, 75014 Paris, France

⁸Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

⁹Department of Astronomy, Columbia University, New York, NY 10027, USA

¹⁰School of Physics and Technology, Wuhan University, Wuhan 430072, People's Republic of China

¹¹Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan

¹²Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People's Republic of China

Received 2024 February 29; revised 2024 December 25; accepted 2025 January 2; published 2025 February 21

## Abstract

JWST is revealing a remarkable new population of high-redshift ($z \gtrsim 4$), low-luminosity active galactic nuclei in deep surveys and detecting the host galaxy's stellar light in the most luminous and massive quasars at $z \sim 6$ for the first time. Recent findings claim that supermassive black holes (SMBHs) in these systems are significantly more massive than predicted by the local black hole (BH) mass-stellar mass ($M_{BH}-M_{\star}$) relation and that this is not due to sample selection effects. Through detailed statistical modeling, we demonstrate that the coupled effects of selection biases (i.e., finite detection limit and requirements for detecting broad lines) and measurement uncertainties can largely explain the reported offset and flattening in the observed $M_{BH}-M_{\star}$ relation toward the upper envelope of the local relation, even for those at $M_{BH} < 10^8 M_{\odot}$. We further investigate the possible evolution of the $M_{BH}-M_{\star}$ relation at $z \gtrsim 4$ with careful treatment of observational biases and consideration of the degeneracy between intrinsic evolution and dispersion in this relation. The bias-corrected intrinsic $M_{BH}-M_{\star}$ relation in the low-mass regime ($M_{\star} \lesssim 10^{10} M_{\odot}$) suggests a large population of low-mass BHs ($M_{BH} \lesssim 10^5 M_{\odot}$), possibly originating from lighter seeds, may remain undetected or unidentified. These results underscore the importance of forward modeling observational biases to better understand BH seeding and SMBH-galaxy coevolution mechanisms in the early universe, even with the deepest JWST surveys.

*Unified Astronomy Thesaurus concepts: Active galactic nuclei (16); Quasars (1319); Supermassive black holes (1663); Galaxy evolution (594); Early universe (435); Active galaxies (17)*

## 1. Introduction

The James Webb Space Telescope (JWST) is revolutionizing the study of high-redshift active galactic nuclei (AGNs) and the early coevolution of supermassive black holes (SMBHs) with their host galaxies. The unprecedented sensitivity and resolution of NIRCam imaging have enabled the first detection of stellar light from AGNs during the reionization epoch ($z \sim 6$) and estimates of their stellar masses. With NIRSpec, black hole (BH) mass estimates, based on the broad H$\alpha$ and H$\beta$ lines, at such high redshifts are routinely measured with ease. Consequently, several recent studies have investigated the $M_{BH}-M_{\star}$ relation at early cosmic times (e.g., X. Ding et al. 2023; Y. Harikane et al. 2023; D. D. Kocevski et al. 2023; R. Maiolino et al. 2023; M. A. Stone et al. 2024; M. Yue et al. 2024). This relation, with its slope, normalization, and intrinsic scatter being sensitive to the merger history of galaxies and the impact of AGN feedback, has emerged as one of the most important scaling relations in unraveling the nature of SMBH-galaxy

coevolution (e.g., C. Y. Peng 2007; M. Volonteri & P. Natarajan 2009; M. Hirschmann et al. 2010; K. Jahnke & A. V. Macciò 2011; J. Kormendy & L. C. Ho 2013; J. Li et al. 2021a; X. Ding et al. 2022; M. Habouzit et al. 2022; J. I. H. Li et al. 2023; M.-Y. Zhuang & L. C. Ho 2023). Additionally, the $M_{BH}-M_{\star}$ relation, close to the seeding epochs ($z \sim 10-30$), may shed light on the formation mechanisms of SMBHs (e.g., M. Volonteri & P. Natarajan 2009; M. Volonteri 2010; A. Trinca et al. 2022; M. T. Scoggins et al. 2023).

In contrast to expectations based on the tight $M_{BH}-M_{\star}$ relation in the local universe (e.g., J. Kormendy & L. C. Ho 2013), AGNs currently found at $4 < z < 7$ tend to significantly lie above this local relation and exhibit larger scatter in the $M_{BH}-M_{\star}$ plane. The detection of “overmassive” BHs (i.e., elevated above the mean local $M_{BH}-M_{\star}$ relation), at least to some extent, is anticipated given a well-understood selection bias (T. R. Lauer et al. 2007) since only luminous and massive BHs with emission sufficiently dominating their host galaxy can be observed in a flux-limited survey (e.g., T. Treu et al. 2007; A. Schulze & L. Wisotzki 2011; M. Volonteri & D. P. Stark 2011; Y. Shen et al. 2015; M. Sun et al. 2015; J. Li et al. 2021a, 2022; M. Volonteri et al. 2023; H. Zhang et al. 2023; T. S. Tanaka et al. 2025). This selection bias is suggested

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
to potentially explain the $\gtrsim 1$ dex offset (relative to the local relation) found for $z \sim 6$ quasars by A. Schulze & L. Wisotzki (2014) and J. Li et al. (2022), although their analyses are based on [C II]-traced dynamical mass measured by the Atacama Large Millimeter/submillimeter Array (ALMA). Contrary to this, several recent JWST studies have argued that selection bias is either negligible in deep JWST surveys or insufficient to explain the large observed offset at high redshift (e.g., F. Pacucci et al. 2023; R. Maiolino et al. 2023; M. A. Stone et al. 2024).

In particular, F. Pacucci et al. (2023) estimate that the observed $M_{BH}-M_*$ relation at $z \sim 5$ for low-luminosity AGNs discovered in the deep JADES (D. J. Eisenstein et al. 2023) and CEERS (S. L. Finkelstein et al. 2023) surveys is a factor of $\sim 10–100$ above that of local AGNs and the bias due to selection is 0.2 dex upward at most. Consequently, they conclude that the overall population of high-redshift BHs is intrinsically overmassive compared to their low-$z$ counterparts. Reaching a similar conclusion, M. A. Stone et al. (2024) compare the distribution of $z \sim 6$ quasars to 30 Palomar-Green (PG) quasars ($z < 0.5$) and find that these two samples occupy distinct regions in the $M_{BH}-M_*$ plane. This is treated as evidence that the elevated BH mass seen in high-redshift quasars cannot be attributed to selection effects.

If selection effects were minimal, particularly at the lower luminosities probed by deep JWST fields, the prevalence of overmassive observed AGNs during the reionization epoch suggests that the majority of the high-redshift BH population is overmassive. This finding could support the theory that SMBHs mainly originate from heavy seeds via direct collapse BHs ($\sim 10^5 M_\odot$) (e.g., M. C. Begelman et al. 2006; B. Agarwal et al. 2013; E. Visbal & Z. Haiman 2018; F. Pacucci et al. 2023; Á. D. Bogdán et al. 2024), with SMBHs initially much more massive than the budding host galaxy (B. Agarwal et al. 2013).

However, the impact of selection effects, measurement uncertainties, and appropriate choice of the local relation is not yet well established to exclude the existence of undermassive BHs (i.e., located on or below the local relation) grown from light seeds ($\sim 100 M_\odot$, formed from the remnants of massive Population III stars) (C. L. Fryer et al. 2001) or medium-weight seeds ($\sim 1000 M_\odot$, formed from runaway stellar collisions in dense star clusters) (P. Madau & M. J. Rees 2001) since they may not have grown enough to be detectable (P. Natarajan et al. 2017; M. Volonteri et al. 2023). Specifically, F. Pacucci et al. (2023) adopted a constant detection limit for $M_{BH}$, derived from the $H\alpha$ line sensitivity, in their bias analysis. However, the detection limit on $M_{BH}$ in a flux-limited survey should depend on the Eddington ratio ($\lambda_{Edd}$) distribution function (ERDF) of the AGN population. A lower intrinsic $\lambda_{Edd}$ would result in only massive BHs being observable. Of similar concern, the PG quasars, used as a low-redshift benchmark in M. A. Stone et al. (2024), have lower luminosities and a different selection function compared to high-redshift quasars, due to being drawn from different mass, accretion rate, and duty cycle distribution functions at low and high redshifts (A. Schulze & L. Wisotzki 2011).

In this paper, we conduct a detailed statistical modeling of selection biases for a sample of high-redshift AGNs observed by JWST, with a focus on recently discovered low-luminosity AGNs, to understand the impact of selection effects on revealing the intrinsic connection between SMBH growth and

galaxy evolution (i.e., the $M_{BH}-M_*$ relation) in the early universe ($z \sim 4–7$). Furthermore, we address the use of different local relations as the baseline, particularly between active and non-active samples. Throughout this paper, we adopt a flat $\Lambda$CDM cosmology with $\Omega_\Lambda = 0.7$ and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. Magnitudes are given in the AB system. Stellar masses are matched to a G. Chabrier (2003) initial mass function (IMF).

## 2. Sample

In this work, we analyze 32 AGNs at $z > 4$, including luminous quasars, from the following samples with JWST observations.

1. Twelve low-luminosity AGNs as reported by R. Maiolino et al. (2023),<fnref n="13" /> and 10 low-luminosity AGNs from Y. Harikane et al. (2023) at $4 < z < 7$. These low-luminosity AGNs, with bolometric luminosities and BH masses in the range of $44.1 \lesssim \log L_{bol}/\text{erg s}^{-1} \lesssim 46.0$ and $6.5 \lesssim \log M_{BH}/M_\odot \lesssim 8.3$, were recently discovered in deep JWST surveys (hereafter the JDEEP sample), including the CEERS (S. L. Finkelstein et al. 2023), GLASS (T. Treu et al. 2022), and JADES (D. J. Eisenstein et al. 2023) programs. They are selected based on the detection of broad components in the $H\alpha$ or $H\beta$ lines from a parent galaxy sample observed with NIRSpec. This selection typically requires a full width at half-maximum (FWHM) $\gtrsim 1000 \text{ km s}^{-1}$, along with additional criteria to ensure the significance of line detection and to exclude cases where line broadening is caused by outflows. To ensure a relatively uniform selection of BL AGNs and a reasonable understanding of the underlying distribution functions for statistical analysis (see Section 4), we excluded individual high-$z$ AGNs reported in the literature or those selected by alternative criteria (e.g., X-ray detection or V-shaped little red dots) with fundamentally different selection functions (e.g., A. D. Goulding et al. 2023; V. Kokorev et al. 2023, 2024; Á. D. Bogdán et al. 2024; J. E. Greene et al. 2024; H. Übler et al. 2024).

2. Two moderate-luminosity quasars with $\log L_{bol}/\text{erg s}^{-1} \sim 46.2$ and $\log M_{BH}/M_\odot \sim 8.3 – 9.1$ from the SHELLQs survey (Y. Matsuoka et al. 2016) as reported in X. Ding et al. (2023). These two quasars are part of the JWST Cycle 1 program (Observation ID 1967; PI: M. Onoue), which targets 12 of the moderate-luminosity quasars with $45.8 \lesssim \log L_{bol}/\text{erg s}^{-1} \lesssim 46.4$ and $7.5 \lesssim \log M_{BH}/M_\odot \lesssim 9.1$ at $6.0 < z < 6.4$ discovered by the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP; H. Aihara et al. 2022).

3. Ten luminous quasars in M. A. Stone et al. (2024) and M. Yue et al. (2024) at $5 < z < 7$, which probe the most luminous and massive BHs in the early universe ($\log L_{bol}/\text{erg s}^{-1} \gtrsim 47.2$ and $\log M_{BH}/M_\odot \gtrsim 8.3$). These quasars are specifically selected to cover a range of luminosity, star formation rate, and the properties of the Ly$\alpha$ forest and Gunn-Peterson trough to study the early growth of SMBHs, intergalactic medium, and

<footnote n="13">¹³ We do not include the secondary BHs of the three dual AGN candidates identified through two broad $H\alpha$ components, as their stellar mass is unknown and their nature is yet to be confirmed.</footnote>
cosmic reionization. Note that one of the quasars in M. A. Stone et al. (2024), HSC J2239+0207, is selected from the SHELLQs survey, with its luminosity (log L<sub>bol</sub>/erg s<sup>-1</sup> ~ 46.4) being an order of magnitude smaller than the rest of the sample. Therefore, we attribute it to the SHELLQs sample in the following analysis.

The BH masses of these AGNs are derived using the single-
epoch virial estimator from the broad Hα or Hβ line, with the
M. Vestergaard & B. M. Peterson (2006) calibration adopted
for Hβ-based BH masses. For consistency, we recalibrated the
Hα-based BH mass in Y. Harikane et al. (2023) using the
updated A. E. Reines & M. Volonteri (2015) recipe, as adopted
in R. Maiolino et al. (2023). This calibration aligns with the
Hβ-based M<sub>BH</sub> from M. Vestergaard & B. M. Peterson (2006)
and the latest reverberation mapping results in Y. Shen et al.
(2024a). In contrast, the J. E. Greene & L. C. Ho (2005) recipe
used in Y. Harikane et al. (2023) yields an Hα-based M<sub>BH</sub> that
is systematically smaller by ~0.3 dex compared to the Hβ-
based calibration.

The L<sub>bol</sub> of low-luminosity AGNs are derived from the broad Hα luminosity (L<sub>Hα</sub>) using the L<sub>Hα</sub>-L<sub>5100</sub> relation in J. E. Greene & L. C. Ho (2005), assuming L<sub>bol</sub> = 10.3 × L<sub>5100</sub> (G. T. Richards et al. 2006; A. E. Reines & M. Volonteri 2015). This approach differs from Y. Harikane et al. (2023), who used the broad+narrow Hα luminosity to derive L<sub>bol</sub>, and from R. Maiolino et al. (2023), who used the L<sub>Hα</sub>-L<sub>bol</sub> relation in J. Stern & A. Laor (2012). For luminous quasars, the contamination of the host galaxy to L<sub>5100</sub> is negligible, allowing their L<sub>bol</sub> to be directly estimated from the continuum luminosity at 5100 Å.

Stellar masses are estimated through various methods: NIRSpec (plus NIRCam photometry) spectral decomposition in R. Maiolino et al. (2023), 2D AGN+host image decomposition with NIRCam data in Y. Harikane et al. (2023), X. Ding et al. (2023), and M. Yue et al. (2024), as well as 1D surface brightness profile scaling in M. A. Stone et al. (2024). Table 1 summarizes the updated masses and luminosities for the JDEEP sample adopted in this work.

We assume that these BH and stellar mass estimates are
unbiased on average, with a statistical uncertainty of ~0.4 dex.
However, in practice, the virial M<sub>BH</sub> may be biased high,
especially for highly accreting systems (Y. Shen 2013; P. Du
et al. 2016). The stellar mass estimates are likely biased due to
systematic effects and degeneracies in spectra- or imaging-
based AGN-host decomposition, especially for AGN-domi-
nated systems at faint magnitudes (e.g., J. Li et al. 2021b;
J. I. H. Li et al. 2023; J. F. W. Baggen et al. 2024;
M.-Y. Zhuang & Y. Shen 2024), imperfect alignment of
objects in the micro-shutter leading to substantial slit loss (e.g.,
S. Tacchella et al. 2023), older stars being overwhelmed by
ongoing star formation (e.g., M. W. Topping et al. 2022;
C. Giménez-Arteaga et al. 2023, 2024; D. Narayanan et al.
2024), and uncertainties in star formation history and IMF
at early cosmic times (e.g., C. L. Steinhardt et al.
2022; M. W. Topping et al. 2022; L. Whitler et al. 2023;
C. Giménez-Arteaga et al. 2024). Such biases are tentatively
indicated by the distinct M<sub>*</sub> distributions for the Y. Harikane
et al. (2023) and R. Maiolino et al. (2023) samples, with a
p-value = 0.007 in a Kolmogorov-Smirnov test and a mean
(median) M<sub>*</sub>difference of 0.43 (0.75) dex, despite having
similar distributions in M<sub>BH</sub> and L<sub>bol</sub>.

<table><caption>**Table 1**
Revised Masses and Luminosities of the JDEEP Sample</caption>
  <thead>
    <tr>
      <th>ID</th>
      <th>z</th>
      <th>M<sub>*</sub><br>(M<sub>☉</sub>)</th>
      <th>M<sub>BH</sub><br>(M<sub>☉</sub>)</th>
      <th>L<sub>bol</sub><br>(erg s<sup>-1</sup>)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>8083</td>
      <td>4.648</td>
      <td>8.45<sup>+0.10</sup><sub>-0.10</sub></td>
      <td>7.28<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>44.81</td>
    </tr>
    <tr>
      <td>1093</td>
      <td>5.595</td>
      <td>8.34<sup>+0.20</sup><sub>-0.20</sub></td>
      <td>7.38<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>44.98</td>
    </tr>
    <tr>
      <td>3608</td>
      <td>5.269</td>
      <td>8.38<sup>+0.15</sup><sub>-0.15</sub></td>
      <td>6.83<sup>+0.40</sup><sub>-0.40</sub></td>
      <td>44.29</td>
    </tr>
    <tr>
      <td>11836</td>
      <td>4.409</td>
      <td>7.79<sup>+0.30</sup><sub>-0.30</sub></td>
      <td>7.11<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>44.72</td>
    </tr>
    <tr>
      <td>20621</td>
      <td>4.681</td>
      <td>8.06<sup>+0.70</sup><sub>-0.70</sub></td>
      <td>7.32<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>44.90</td>
    </tr>
    <tr>
      <td>77652</td>
      <td>5.229</td>
      <td>7.87<sup>+0.28</sup><sub>-0.28</sub></td>
      <td>6.84<sup>+0.39</sup><sub>-0.39</sub></td>
      <td>44.72</td>
    </tr>
    <tr>
      <td>61888</td>
      <td>5.875</td>
      <td>8.11<sup>+0.92</sup><sub>-0.92</sub></td>
      <td>7.21<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>44.98</td>
    </tr>
    <tr>
      <td>62309</td>
      <td>5.172</td>
      <td>8.12<sup>+0.13</sup><sub>-0.13</sub></td>
      <td>6.54<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>44.46</td>
    </tr>
    <tr>
      <td>954</td>
      <td>6.760</td>
      <td>10.66<sup>+0.10</sup><sub>-0.10</sub></td>
      <td>7.89<sup>+0.35</sup><sub>-0.35</sub></td>
      <td>45.67</td>
    </tr>
    <tr>
      <td>10013704</td>
      <td>5.919</td>
      <td>8.88<sup>+0.66</sup><sub>-0.66</sub></td>
      <td>7.48<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>44.55</td>
    </tr>
    <tr>
      <td>73488</td>
      <td>4.133</td>
      <td>9.78<sup>+0.20</sup><sub>-0.20</sub></td>
      <td>7.71<sup>+0.35</sup><sub>-0.35</sub></td>
      <td>45.16</td>
    </tr>
    <tr>
      <td>53757</td>
      <td>4.448</td>
      <td>10.18<sup>+0.13</sup><sub>-0.13</sub></td>
      <td>7.52<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>44.64</td>
    </tr>
    <tr>
      <td>CEERS_01244</td>
      <td>4.478</td>
      <td>8.63<sup>+1.03</sup><sub>-1.03</sub></td>
      <td>7.73<sup>+0.35</sup><sub>-0.35</sub></td>
      <td>45.16</td>
    </tr>
    <tr>
      <td>GLASS_160133</td>
      <td>4.015</td>
      <td>&lt;8.82</td>
      <td>6.66<sup>+0.35</sup><sub>-0.35</sub></td>
      <td>44.46</td>
    </tr>
    <tr>
      <td>GLASS_150029</td>
      <td>4.583</td>
      <td>9.10<sup>+0.37</sup><sub>-0.37</sub></td>
      <td>6.88<sup>+0.35</sup><sub>-0.35</sub></td>
      <td>44.32</td>
    </tr>
    <tr>
      <td>CEERS_00746</td>
      <td>5.624</td>
      <td>&lt;9.11</td>
      <td>7.91<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>45.97</td>
    </tr>
    <tr>
      <td>CEERS_01665</td>
      <td>4.483</td>
      <td>9.92<sup>+0.68</sup><sub>-0.68</sub></td>
      <td>7.51<sup>+0.35</sup><sub>-0.35</sub></td>
      <td>45.11</td>
    </tr>
    <tr>
      <td>CEERS_00672</td>
      <td>5.666</td>
      <td>&lt;9.01</td>
      <td>7.90<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>45.48</td>
    </tr>
    <tr>
      <td>CEERS_02782</td>
      <td>5.241</td>
      <td>&lt;9.35</td>
      <td>7.85<sup>+0.36</sup><sub>-0.36</sub></td>
      <td>45.15</td>
    </tr>
    <tr>
      <td>CEERS_00397</td>
      <td>6.000</td>
      <td>9.36<sup>+0.45</sup><sub>-0.45</sub></td>
      <td>7.29<sup>+0.35</sup><sub>-0.35</sub></td>
      <td>44.75</td>
    </tr>
    <tr>
      <td>CEERS_00717</td>
      <td>6.936</td>
      <td>9.61<sup>+0.18</sup><sub>-0.18</sub></td>
      <td>8.28<sup>+0.40</sup><sub>-0.40</sub></td>
      <td>44.46</td>
    </tr>
    <tr>
      <td>CEERS_01236</td>
      <td>4.484</td>
      <td>8.94<sup>+0.54</sup><sub>-0.54</sub></td>
      <td>7.58<sup>+0.35</sup><sub>-0.35</sub></td>
      <td>44.11</td>
    </tr>
  </tbody>
</table>

Note. The BH mass uncertainty is derived by adding 0.35 dex, which represents the statistical uncertainty of the virial estimator, in quadrature to the measurement uncertainties from emission line width and luminosity. We impose a minimum stellar mass uncertainty of 0.1 dex for values below this threshold. The maximum value of the asymmetric error bar of each object is adopted in this work following F. Pacucci et al. (2023).

It is also noteworthy that, although the stellar masses for
many “undetected” host galaxies are reported as upper limits,
these systems might actually be “massive” (compared to the
upper limits) but compact galaxies indistinguishable from a
point source (e.g., J. Li et al. 2021b; J. F. W. Baggen et al.
2024; M.-Y. Zhuang & Y. Shen 2024). This effect is
particularly pronounced for quasar-dominated systems and
barely resolved low-mass (M<sub>*</sub> < 10<sup>9</sup> M<sub>☉</sub>) galaxies at z ≳ 5,
where the typical galaxy size for the latter (≈0.5 kpc, ≲0.09) is
smaller than the PSF in most long wavelength filters
(J. F. W. Baggen et al. 2023; R. Maiolino et al. 2023).
Dedicated image and spectral energy distribution simulations
are needed to understand how reliably the host properties can
be recovered in such systems (J. Li et al. 2021b; M.-Y. Zhuang
& Y. Shen 2024).

With these caveats in mind, we present the M<sub>BH</sub>−M<sub>*</sub>
relation for our compiled high-redshift AGN sample in
Figure 1. The distinct distributions of various samples
demonstrate how surveys with varying areas and depths
complicate our understanding of this relation. The rare, high-
luminosity quasars in massive host galaxies come from
individual pointings preselected from all-sky but shallow
surveys (X. Fan et al. 2023), which are only capable of
discovering the most massive and luminous systems. The wide
and deep HSC-SSP survey has extended the detection limit to
somewhat lower BH masses (Y. Matsuoka et al. 2016;
M. Onoue et al. 2019), while the very low-luminosity AGNs
<figure><img src="image_5.png" /><figcaption>**Figure 1.** The $M_{BH}-M_{\star}$ relation of high-redshift AGNs used in this study. The filled stars, diamonds, squares, circles, and triangles represent AGNs in R. Maiolino et al. (2023), Y. Harikane et al. (2023), X. Ding et al. (2023), M. A. Stone et al. (2024), and M. Yue et al. (2024), respectively, with each point color coded by $L_{bol}$. For comparison, we show the SDSS quasar sample ($\log L_{bol}/\text{erg s}^{-1} \lesssim 46$) at $0.2 < z < 1$ from J. Li et al. (2021a, 2021b) as blue contours (1$\sigma$-3$\sigma$ level), along with the local $M_{BH}-M_{\text{bulge}}$ relation from J. Kormendy & L. C. Ho (2013) and the local $M_{BH}-M_{\star}$ relation from J. E. Greene et al. (2020). The gray shaded regions represent a schematic view of the inaccessible regions of high-redshift observations, where BHs are either undetectable or underrepresented (sparsely sampled) due to limited sensitivity and survey volume, with the exact limits depending on specific survey designs.</figcaption></figure>

are detected in deep but pencil-beam JWST fields that mainly probe common low-mass galaxies. This explains why, at $\log M_{BH}/M_{\odot} > 9.5$, the AGNs identified in deep JWST surveys exhibit significantly lower luminosities and $M_{BH}/M_{\star}$ ratios (by $\sim 1.5$ dex) compared to those in quasar samples. At the high-$M_{\star}$ end, these AGNs are the most common objects in deep pencil-beam surveys, while quasar-like objects, although detectable, are much rarer.

### 3. $M_{BH}-M_{\star}$ Relation at High-$z$

To elucidate the connection between the early growth of SMBHs and their host galaxies, a common approach is to compare their positions in the $M_{BH}-M_{\star}$ plane with those observed in the local universe. In this study, we adopt the following local relations (Figure 1): (1) The $M_{BH}-M_{\text{bulge}}$ relation with an intrinsic scatter of $\sim 0.3$ dex from J. Kormendy & L. C. Ho (2013, hereafter KH13). This relation, established with massive (classical) bulge-dominated quiescent galaxies with dynamical BH mass measurements, has been a standard for exploring SMBH–galaxy coevolution. Note that KH13 increased the mass-to-light ratio by 0.125 dex to match the dynamical zero-point when measuring the bulge mass. We use their unadjusted relation to maintain consistency with our AGN sample. (2) The $M_{BH}-M_{\star}$ relation for “all” local massive BHs ($2 < \log M_{BH} < 10$) from J. E. Greene et al. (2020, hereafter G20), established through a combined fit of local early- and late-type galaxies with (mostly) dynamical BH mass measurements. This relation, based on total stellar mass rather than bulge mass, exhibits a much lower normalization and larger scatter ($\sim 0.81$ dex) compared to the KH13 relation.

While the observed high-redshift AGNs do exhibit an elevated $M_{BH}/M_{\star}$ compared to the mean local ratios shown in Figure 1, the dependence of $M_{BH}$ on $L_{bol}$ ensures that the

T. R. Lauer et al. (2007) bias, induced by flux limits and the intrinsic scatter of the mass relation must be present. There will always be regions in the high-$z$ $M_{BH}-M_{\star}$ plane that are either not accessible due to the finite detection limits, or sparsely populated due to the combined effects of the limited survey area, the scatter of the mass relation, and the shape of the stellar mass function (SMF). A schematic view of these inaccessible regions is displayed in Figure 1 as gray shaded regions. The exact limits depend on the specific sample and can vary strongly with population and redshift.

In Figure 1, we also show the $M_{BH}-M_{\star}$ relation for $\sim 4000$ Sloan Digital Sky Survey (SDSS) quasars at $0.2 < z < 1.0$ and above a redshift-dependent stellar mass limit from J. Li et al. (2021a, 2021b) as blue contours. In contrast to the small sample comparison in M. A. Stone et al. (2024), we find that the distribution of high-redshift quasars is covered by the upper envelope of the SDSS quasar sample. Therefore, “over-massive” BHs exist in both high- and low-redshift universe. This is true even considering that the low-$z$ quasar sample in J. Li et al. (2021b) probes only $\log L_{bol}/\text{erg s}^{-1} \lesssim 46$ and misses the most luminous and massive quasars that are sparse on the sky, due to their requirement of having HSC-SSP DR2 coverage ($\sim 300 \text{ deg}^2$) for measuring stellar masses. The absence of $z \sim 6$ quasars populating the lower envelope of the SDSS sample may be due to the shallower detection limit at high redshift, resulting in the undetection of the entire “undermassive” population. In the following sections, we conduct a rigorous modeling of observational biases to support this hypothesis.

### 4. Method

We quantify the impact of selection effects and measurement uncertainties on the intrinsic mass relation by modeling the conditional probability of observing a $M_{BH}$ for a given $M_{\star}$ at a specified detection limit. We then infer the intrinsic mass relation from the observed biased sample through forward modeling as widely adopted in the literature (e.g., A. Schulze & L. Wisotzki 2011; Y. Shen et al. 2015; M. Sun et al. 2015; J. Li et al. 2021a, 2022; X. Ding et al. 2023; J. I. H. Li et al. 2023; T. S. Tanaka et al. 2025). Here, we outline the statistical recipes for measuring the bivariate distribution of the observed $M_{BH}$ and $M_{\star}$.

We define the true mass terms as $m \equiv \log M_{BH}$ and $s \equiv \log M_{\star}$, while the observed mass measurements, including Gaussian uncertainties $\sigma_m$ and $\sigma_s$, are denoted as $m_o$ and $s_o$, respectively. The luminosity of the quasar $l \equiv \log L_{bol}$ is related to $m$ through the Eddington ratio $\lambda \equiv \log \lambda_{Edd}$ as $l = \lambda + m + 38.1$. We ignore the uncertainty on $l$, which is on the order of $\sim 0.04$ dex (Y. Harikane et al. 2023; R. Maiolino et al. 2023); although we note that this uncertainty only accounts for errors propagated from the emission line width and continuum flux, and not the bolometric conversion factor. The observed Eddington ratio $\lambda o$ can then be measured as $l - m_o - 38.1$.

At a given true $m$, the line width ($w \equiv \log \text{FWHM}$, assuming negligible uncertainty) of the broad line (BL) follows a Gaussian distribution $g(w)$ with the mean determined by the virial estimator (M. Vestergaard & B. M. Peterson 2006) as

$$w = 0.5 m - 0.25 l + 10.8. \quad (1)$$

Its scatter $\sigma_w$ (0.14 dex), along with the scatter in $l$ at a fixed $m$ (i.e., the scatter of the ERDF, which is about 0.4 dex;
J. Wu et al. 2022), contributes to the uncertainty in the virial BH mass estimate (0.35 dex) as

$$ \sigma_m = \sqrt{(0.5\sigma_l)^2 + (2\sigma_w)^2}. \qquad (2) $$

We assume that *m* is correlated with *s* through a linear relation *m* = *αs* + *β* with a slope *α*, normalization *β*, and a mass-independent Gaussian intrinsic scatter *σ*, i.e.,

$$ g(m|s) = \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{(m - (\alpha s + \beta))^2}{2\sigma^2}\right). \qquad (3) $$

Note that for mass relations with a steep slope, small normalization, or large scatter, the BH mass below the mean relation for low-mass galaxies could become negative. In such cases, the distribution of $M_{BH}$ at a given $M_*$ can be equivalently considered as a truncated Gaussian distribution, where values of $M_{BH} < 100 M_\odot$ are set to $M_{BH} = 100 M_\odot$. This mass distribution could represent the presence of a significant population of light seeds that fail to grow efficiently (e.g., B. D. Smith et al. 2018).

The probability of detecting an AGN at a redshift of *z*, with true BH mass *m*, stellar mass *s*, and luminosity *l*, can be expressed as

$$ p(m, s, l, z) = \psi(s, z) \times \frac{dV}{dz} \times g(m|s, z) \times p_{ac}(m, z) \times p(\lambda, z) \times \Omega(m, \lambda, w, z), \qquad (4) $$

where $\psi(s, z)$ is the SMF, $p_{ac}$ is the active unobscured BL AGN fraction, $p(\lambda, z)$ is the ERDF, and $\Omega(m, \lambda, w, z)$ is the selection function of observations.

We attribute the selection function primarily to the AGN luminosity limit ($l_{lim}$) of a specific survey and the requirement of having line widths $\geq 1000$ km s$^{-1}$ (e.g., Y. Harikane et al. 2023), which inevitably excludes low-mass BHs with intrinsically narrower BLs. The luminosity limit is determined as the maximum value between $L_{phot,lim}$ and $L_{spec,lim}$. Here, $L_{phot,lim}$ is the minimum $L_{bol}$ above which an AGN can be photometrically selected for scheduling spectroscopic follow-up observations (e.g., the SHELLQs and quasar samples). The $L_{spec,lim}$ is determined by the emission line sensitivity of spectroscopic observations and the requirement for robustly identifying BLs, which can be converted into a detection limit on $L_{bol}$ via its tight correlation with line luminosity (J. E. Greene & L. C. Ho 2005).

For simplicity, we use the lowest luminosity observed for each sample ($\log L_{bol}/\text{erg s}^{-1} = 44.1, 46.2$, and $47.2$ for the JDEEP, SHELLQs, and quasar samples, respectively) as the representative luminosity limit. A more detailed treatment of the selection function would account for host galaxy contamination, complex photometric selection criteria, target prioritization in shutter allocation (which may bias toward or against AGNs) (R. Maiolino et al. 2023), and additional requirements for line detections and specific science goals. However, such an analysis is beyond the scope of the current work. Our approach of adopting an effective luminosity limit based on that of spectroscopically confirmed AGNs effectively reflects the current capability to detect and identify BL AGNs with specific survey designs and color selection methods.

In our model, BHs can be detected as BL AGNs if their luminosity *l* exceeds the detection limit $l_{lim}$ and the line width

exceeds $1000$ km s$^{-1}$. Therefore, the selection function on *m* can be expressed as

$$ \Omega(m) = p_{ac}(m) g_w(w > 3|l| > l_{lim}) \times \int_{l_{lim}}^{l_{max}} p_{\lambda}(l - m - 38.1) dl. \qquad (5) $$

Given the lack of observational constraints, we assume that $p_{ac}$ does not vary with $M_{BH}$, and thus it can be removed from the above Equation. The potential impact of a mass-dependent $p_{ac}$ is discussed in the Appendix. In general, a positive correlation between $p_{ac}$ and $M_{BH}$ at high redshift (e.g., due to increased fueling efficiency) tends to bias the observed mass relation of AGNs toward massive BHs, even without any selection on luminosity. This contrasts with the situation in the local universe, where less massive BHs are preferentially observed in AGN samples due to reduced accretion activity in the most massive BHs (A. Schulze & L. Wisotzki 2011).

The bivariate distribution of the observed $m_o$ and $s_o$, incorporating measurement uncertainties and the selection function, is formulated as

$$ p(m_o, s_o) = \int g(m_o|m) g(s_o|s) p(m, s, z) dm ds. \qquad (6) $$

The observed $M_{BH}-M_*$ relation (i.e., the expectation value of $m_o$ at a given $s_o$) can then be derived as

$$ m_o(s_o) = \frac{\int m_o p(m_o, s_o) dm_o}{\int p(m_o, s_o) dm_o}. \qquad (7) $$

Measuring the bivariate $M_{BH}-M_*$ distribution and the expected “observed” mass relation given the selection function and mass uncertainties relies on prior knowledge of the underlying intrinsic distribution functions (see Equations (4), (6), and (7)). In this work, we adopt the latest SMF from R. Navarro-Carrera et al. (2024), which utilized deep JWST data to constrain the slope at the low-mass end. For the ERDF, we use the fiducial result from J. Wu et al. (2022), which can be approximated by a log-normal distribution with a mean of $\sim -1.0$ and a scatter of $\sim 0.4$. This ERDF was derived for quasars at $z \sim 6$ with selection effects corrected through a Bayesian framework.

The BL AGN fraction $p_{ac}$ remains uncertain at $z \sim 5$. In the $H\alpha$ flux-limited EIGER and FRESCO NIRCam/WFSS surveys, the fraction of BL AGNs with $L_{H\alpha} > 2.3 \times 10^{42} \text{erg s}^{-1}$ ($\log L_{bol}/\text{erg s}^{-1} > 44.7$) to UV-selected galaxies is $\sim 1\%$ at $-19 < M_{UV} < -22$ ($8 \lesssim \log M_*/M_\odot \lesssim 10$; M. Song et al. 2016). This fraction increases to $\sim 5\%$ in the deeper CEERS + GLASS surveys (Y. Harikane et al. 2023), and further rises to $\sim 10\%$ in the deepest JADES survey (R. Maiolino et al. 2023), although it may be affected by the complex survey selection function. Conservatively, the fraction of all unobscured BL AGNs considered by our model (i.e., including those below the flux limit) must exceed $\sim 1\%$. On the other hand, the AGN duty cycle at $z \sim 6$, deduced from solving the continuity equation, ranges from $\sim 15\%$ at $\log M_{BH}/M_\odot \sim 6.5$ and increases to $60\%$ for $\log M_{BH}/M_\odot \sim 8.3$ (R. Aversa et al. 2015). This duty cycle can be converted into the unobscured BL AGN fraction by assuming a ratio of unobscured to obscured AGNs. Despite substantial uncertainties and sample incompleteness, this ratio is $\sim 1:2$ in the JADES survey (R. Maiolino et al. 2023; J. Scholtz et al. 2023), which is broadly consistent with constraints from X-ray observations (F. Vito et al. 2018). Based on these findings, we
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Impact of selection effects and measurement uncertainties on the $M_{BH}-M_{\star}$ relation. The square data points, color coded by the average $L_{bol}$ in each mass bin, represent mock AGNs simulated using Equation (6) based on the local KH13 relation (red dashed line) with an intrinsic scatter of 0.3 dex. The cyan solid curve represents the mean $M_{BH}-M_{\star}$ relation averaged over the square data points. Panel (a): The original KH13 relation. Panel (b): no mass uncertainties are included, but the selection function (log $L_{bol}$ > 44.1 erg s⁻¹ and FWHM > 1000 km s⁻¹) is applied to the input KH13 relation. The effects of applying the $L_{bol}$ and FWHM cuts separately are also shown for comparison. Panel (c): mass uncertainties are included for both $M_{BH}$ (0.35 dex) and $M_{\star}$ (0.45 dex), but no selection function is applied. Panel (d): mass uncertainties are included and the selection function is applied.</figcaption></figure>

consider $p_{ac} \sim 1\%–5\%$ for BL AGNs (assumed to be mass independent) in this study. The impact of the assumed $p_{ac}$ range is discussed in Section 5.3.

# 5. Results

## 5.1. Impact of Observational Biases

We first demonstrate how the observational biases (i.e., selection effects induced by flux limits and measurement uncertainties) modify the intrinsic $M_{BH}-M_{\star}$ relation using a Monte Carlo simulation that mimics the effective JDEEP observations. We assume that SMBHs at $z=5.2$ (the mean redshift of the JDEEP sample) follow the local KH13 relation with an intrinsic scatter of 0.3 dex. By randomly sampling Equation (6), we create mock AGNs with known intrinsic $M_{BH}$, $M_{\star}$, $L_{bol}$, and line width. The distribution of the resulting mock sample in the $M_{BH}-M_{\star}$ plane, color coded by

the average $L_{bol}$ in each mass bin, is displayed in Figure 2(a). It is evident that $M_{BH}$ strongly depends on $L_{bol}$, which is the fundamental cause of the selection bias.

In Figure 2(b), we show the mass relation after applying the selection criteria of $\log L_{bol}/\text{erg s}^{-1} > 44.1$ and $\text{FWHM} > 1000 \text{ km s}^{-1}$ to the mock AGNs from Figure 2(a). It exhibits a significant positive offset at the low-$M_{\star}$ end due to missing low-mass BHs. The luminosity cut is primarily responsible for this offset, while the FWHM cut further excludes a portion of bright AGNs with a line width narrower than $1000 \text{ km s}^{-1}$ (Figure 2(b)).

In Figure 2(c), we incorporate the average mass uncertainties for the JDEEP sample ($\sigma_m = 0.35$ and $\sigma_s = 0.45$) into the mock AGNs from Figure 2(a) to illustrate the impact of measurement uncertainties. Notably, the resulting $M_{BH}-M_{\star}$ distribution is significantly broader than the original KH13 relation, and its mean appears to shift downward. This systematic shift, known
<figure><img src="image_7.png" /><figcaption>**Figure 3.** The predicted “observed” $M_{BH}-M_{\star}$ relations (solid curves) for detectable AGNs at $z \sim 5$ in the three surveys, assuming that the underlying BH population intrinsically follows the local KH13 and G20 relations (dashed curves with the same colors). The intrinsic scatter at $1\sigma$, $2\sigma$, and $3\sigma$ levels for each local relation is shown as shaded regions. Only BHs in the upper envelope of the assumed intrinsic mass relations can be detected as luminous AGNs. Their positions are further biased upward by the substantial mass uncertainties (Figure 2), resulting in a significantly elevated $M_{BH}/M_{\star}$ compared to the intrinsic ratio.</figcaption></figure>

as the A. S. Eddington (1913) bias, is caused by the substantial stellar mass uncertainties and the bottom-heavy SMF. To elaborate, consider a stellar mass bin at log $M_{\star}/M_{\odot} \sim 8.5$. Due to the abundance of small BHs hosted by low-mass galaxies with log $M_{\star}/M_{\odot} < 8.5$, the added uncertainties on the stellar mass cause numerous low-mass BHs to scatter rightward into the log $M_{\star}/M_{\odot} \sim 8.5$ bin. This effect outweighs the leftward scattering of the rarer high-mass BHs with log $M_{\star}/M_{\odot} > 8.5$ into the same stellar mass bin. Consequently, the mean observed BH mass at log $M_{\star}/M_{\odot} \sim 8.5$ is biased low.

Finally, we apply the log $L_{bol}/\text{erg s}^{-1} > 44.1$ and FWHM $> 1000 \text{ km s}^{-1}$ criteria to the mock AGNs in Figure 2(c). This gives rise to the observed $M_{BH}-M_{\star}$ relation that incorporates both realistic mass uncertainties and the selection function. As shown in Figure 2(d), it appears noticeably flatter than the KH13 relation, with a significant positive offset at the low-$M_{\star}$ end and a negative offset at the high-$M_{\star}$ end, due to the combined effects of the Lauer bias and the Eddington bias. The comparison between Figures 2(b) and 2(d) demonstrates that substantial mass uncertainties further amplify the observed offset in $M_{BH}/M_{\star}$ at the low-$M_{\star}$ end (up to $\sim 1.5$ dex), rather than simply broadening the distribution.

Our analyses underscore the necessity of accounting for the overall transformation in the observed $M_{BH}-M_{\star}$ distribution (i.e., deviation, broadening, and flattening relative to the intrinsic mass relation) caused by selection effects and measurement uncertainties, even in deep JWST surveys like JADES and CEERS.

## 5.2. JWST-observed Overmassive AGNs

Similar to the analyses in Figure 2, we apply the selection function and the average measurement uncertainties of each observed sample to the local KH13 and G20 relations. This allows us to predict the $M_{BH}-M_{\star}$ relation expected at the effective detection limit of each survey using Equation (7), assuming SMBHs at $z \sim 5$ intrinsically follow the local relations. As shown in Figure 3, the predicted $M_{BH}-M_{\star}$

relations above the luminosity limit deviate significantly from the local relations and generally align with the data points within the observed stellar mass ranges, confirming earlier results based on ALMA (A. Schulze & L. Wisotzki 2014; J. Li et al. 2022). Notably, although the G20 relation has a lower average $M_{BH}/M_{\star}$ compared to the KH13 relation, its predicted observed $M_{BH}/M_{\star}$ could be even higher due to its larger scatter, leading to more observable overmassive BHs. Note that the scatter of the G20 relation might be mass dependent and decrease at the massive end, although this was not considered in their fitting.

These results suggest that even if the underlying $z \sim 5$ BH population follows the local relation on average, we would only observe its upper envelope, namely, overmassive AGNs, due to the finite detection limit of JWST or specific sample selection criteria. Our analysis highlights the importance of properly accounting for the dispersion of the mass relation when studying its evolution, rather than simply comparing averages across cosmic times. Figure 3 also exemplifies that one should avoid plotting a data set against an intrinsic correlation without transforming the latter to match the data set's selection and measurement uncertainties. Failing to do so can easily lead to misinterpretations.

It should be noted that the predicted $M_{BH}-M_{\star}$ relation only represents an idealized scenario for an infinitely large survey, as the survey volume, BL AGN fraction, and normalization of the SMF in Equation (4) cancel out in Equation (7). In practice, only BHs within the survey's accessible stellar mass range (determined by its area and depth) can be observed. Consequently, specific parameter spaces in $[\alpha, \beta, \sigma, p_{ac}, M_{\star}, M_{BH}]$ are unobservable in small-area surveys, meaning the expected number of AGNs above the detection limit may be less than one. In other words, the observed number of AGNs in a survey imposes additional constraints on the underlying intrinsic $M_{BH}-M_{\star}$ relation and other distribution functions. While fully utilizing this constraint in modeling the underlying relations is beyond the scope of this work, we use this argument as a sanity check on the validity of the proposition
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Model-predicted permitted region ($15 < N_{\text{obs}} < 65$ for AGNs with log $L_{\text{bol}} > 44.1$, $7.8 < \log M_{\star} < 10.6$, and $4 < z < 7$) in the $\beta-\sigma$ plane for $\alpha = 1.17$ (KH13) and $\alpha = 1.61$ (G20), respectively. The black dashed and solid curves represent $N = 15$ and $N = 65$, respectively. The black crosses mark the $\beta$ and $\sigma$ values for the KH13 and G20 relations. Each bin is color coded by its average $N_{\text{obs}}$.</figcaption></figure>

that high-redshift SMBHs intrinsically follow the KH13 or G20 relations.

We first use the JDEEP survey as an example to illustrate this point. Specifically, by assuming different forms of the intrinsic $M_{\text{BH}}-M_{\star}$ relation and BL AGN fraction, we count the number ($N_{\text{obs}}$) of observable AGNs (i.e., $\log L_{\text{bol}}/\text{erg s}^{-1} > 44.1$ and $\text{FWHM} > 1000 \text{ km s}^{-1}$) in a survey area of $\sim 100 \text{ arcmin}^2$ (8 NIRSpec pointings for CEERS and 8 pointings for JADES; D. J. Eisenstein et al. 2023; Y. Harikane et al. 2023; R. Maiolino et al. 2023)<fnref n="14" /> within the observed stellar mass range ($7.8 < \log M_{\star}/M_{\odot} < 10.6$) at $4 < z < 7$. If the model produces a number of observable BL AGNs commensurate to observations (i.e., 22 sources in JDEEP), we consider it a permitted model. Otherwise, it is deemed forbidden, meaning that specific combinations of $\alpha, \beta, \sigma$, and $p_{\text{ac}}$ either produce too many bright AGNs that are not observed or fail to produce any observable ones. In practice, we adopt $15 < N_{\text{obs}} < 65$ to identify permitted models. This range is somewhat arbitrary and is determined by first adjusting the actual observed AGN counts by a factor of 1.5 (i.e., $N_{\text{obs}} = 15-33$) to accommodate uncertainties in the galaxy number counts from our adopted SMF and in converting $M_{\text{BH}}$ to $L_{\text{bol}}$ via the assumed ERDF (a higher $\lambda_{\text{Edd}}$ would result in more observable AGNs). We further extend the upper limit by an additional factor of two (i.e., $N_{\text{obs}} = 15-65$) to account for potential missing BL AGNs that exceed the luminosity limit but remain unidentified due to various reasons (e.g., outshined by the host galaxy, restricted shutter arrangement).

As shown in Figure 4, the KH13 relation falls well within the permitted region ($N_{\text{obs}} \sim 37$) with $p_{\text{ac}} = 10\%$. The G20 relation shows a slight deviation ($N_{\text{obs}} \sim 8$) with $p_{\text{ac}} = 10\%$, but a slightly higher (or broader) mass relation or $p_{\text{ac}}$ could bring it back into agreement. Similarly, we count the number of observable quasars with $\log L_{\text{bol}}/\text{erg s}^{-1} > 47.2$ and $9.8 < \log M_{\star}/M_{\odot} < 11.5$ over half of the sky (X. Fan et al. 2023).

across $5 < z < 7$, assuming $p_{\text{ac}} = 10\%$. The predicted number based on the KH13 relation is $N_{\text{obs}} \sim 71$, which aligns well with the few tens of quasars of comparable luminosities discovered thus far (X. Fan et al. 2023). The G20 relation, with a similar mean to the KH13 relation at the massive end but a much larger scatter (0.81 dex), predicts $N_{\text{obs}} \sim 910$. This indicates that the intrinsic scatter for massive galaxies cannot be that large if the norm and slope are those of G20, unless a significant population of unobscured quasars has been missed by current color selection methods.

In general, both the KH13 and G20 relations can generate a sufficient number of observable overmassive AGNs under reasonable assumptions and reproduce the observed flat $M_{\text{BH}}-M_{\star}$ relation given the sample selection function (Figure 3), even if their mean $M_{\text{BH}}/M_{\star}$ are significantly smaller than the observed ratio. This analysis demonstrates that the reported offset in the $M_{\text{BH}}-M_{\star}$ plane can be replicated by incorporating observational biases (including measurement uncertainties) into the local relations (i.e., by sampling their $\gtrsim 2\sigma$ upper envelope), even assuming that the mass estimates are unbiased on average.

### 5.3. Intrinsic $M_{\text{BH}}-M_{\star}$ Relation for JWST-discovered Low-luminosity AGNs

To differentiate whether the observed $M_{\text{BH}}-M_{\star}$ relation is drawn from the KH13, G20, or other forms of intrinsic mass relations that may appear similarly flat after applying the selection function (Figure 3), we employ *emcee* to sample the conditional probability distribution $p(m_o|s_o)$ (i.e., the likelihood function assuming flat priors) derived from Equation (6) as functions of $\theta = [\alpha, \beta, \sigma, p_{\text{ac}}]$. This enables us to infer the most likely intrinsic, bias-corrected $M_{\text{BH}}-M_{\star}$ relation. That is, once convolved with measurement uncertainties and the sample selection function, this relation gives rise to the observed AGN sample.

In this work, we focus on deriving the intrinsic $M_{\text{BH}}-M_{\star}$ relation for the JDEEP sample in the low-$M_{\star}$ regime and comparing our results with those in F. Pacucci et al. (2023). We

<footnote n="14">¹⁴ The effective area for each NIRSpec pointing used for shutter target placement in JADES is ~7 arcmin². Considering pointing overlaps, we use a total area of ~50 arcmin² for JADES and assume the same for CEERS.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Left: best-fit intrinsic $M_{BH}-M_*$ relation (black line) and its $1\sigma$ uncertainty (green shaded region) inferred from the JDEEP sample with solid host galaxy detections (blue points). The inset panel shows the expanded $1\sigma-3\sigma$ intrinsic scatter (gray shaded regions). The intrinsic $M_{BH}-M_*$ relation for the same sample derived by F. Pacucci et al. (2023) and its $1\sigma$ uncertainty are shown in magenta. Local relations from KH13, G20, and RV15 are also shown for comparison. Right: posterior parameter distributions and best-fit values of the mass relation and BL AGN fraction for the observed AGN sample are shown on the left.</figcaption></figure>

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

do not perform a joint fit with the quasar sample, considering that the mass relation may not be linear in the low-$M_*$ and high-$M_*$ regimes (M. Habouzit et al. 2022), and due to the very different selection functions of the two samples. The analysis of the quasar sample will be presented in J. Silverman et al. 2025, in preparation.

In the fitting, we consider only the 18 objects with robust host galaxy detections, as the stellar masses for “undetected” hosts are highly uncertain and not necessarily upper limits (Section 2). In principle, the fitting methodology should self-consistently reproduce the AGN luminosity function given the ERDF, the best-fit $M_{BH}-M_*$ relation, and AGN fraction. However, in practice, this is infeasible due to the small sample size, non-uniform selection function, potential missing BL AGNs, and uncertain underlying distribution functions with unknown mass dependencies. Therefore, we adopt a simplified approach by using the integrated AGN number counts as a constraint in the fitting to eliminate forbidden parameters. Specifically, we define $P(I/\theta) = 1$ if the model predicts $15 < N_{obs} < 65$, and $P(I/\theta) = 0$ otherwise. We use flat priors for all parameters ($-12 < \alpha < 6$, $0.1 < \beta < 2.0$, $0.1 < \sigma < 2.0$, and $1\% < p_{ac} < 15\%$). A selection function of $log L_{bol} > 44.1$ erg s$^{-1}$ and FWHM > 1000 km s$^{-1}$, along with the mass uncertainties for each AGN, were incorporated in the fitting. We then minimize the product of $p(m_o|s_o)$, the prior, and $P(I/\theta)$ as the likelihood function to compute the posterior parameter distributions.

Figure 5 presents the best-fit intrinsic $M_{BH}-M_*$ relation and the corresponding posterior parameter distributions. The best-fit relation is parameterized as $log M_{BH} = 1.02_{-0.54}^{+0.56} \log M_* - 3.32_{-5.16}^{+4.56}$ with an intrinsic scatter of $0.97_{-0.37}^{+0.52}$ dex, corresponding to an observable BL AGN number of $29.3_{-10.8}^{+19.8}$. As expected, the slope and normalization of the intrinsic relation are notably steeper and smaller than the

observed one, after correcting for the selection bias of missing low-mass BHs.

As shown in the inset of Figure 5, the observed AGNs with $log M_*/M_\odot \gtrsim 9.0$ and $log M_{BH}/M_\odot \sim 7.5$ are positioned within the $1\sigma$ intrinsic scatter region. While AGNs with even higher BH masses are detectable, their number density is low, resulting in their underrepresentation in pencil-beam surveys. Conversely, AGNs detected at $log M_*/M_\odot < 9.0$ (where the intrinsic $M_{BH}/M_*$ is lower) lie outside the $1\sigma$ intrinsic scatter region, as only the upper BH envelope is detectable. This result suggests that the overmassive AGNs detected by JWST are consistent with being drawn from an underlying BH population whose average BH mass is much smaller (by $\sim 1.2$ dex) than what has been observed.

The large scatter of $0.97_{-0.37}^{+0.52}$ dex in the low-mass regime, compared to that of the KH13 relation (with a probability of $p(\sigma > 0.3) \sim 98.2\%$) and massive quasars at $z < 2$ ($\sim 0.3$ dex; J. Li et al. 2021a; T. S. Tanaka et al. 2025), may suggest a merger-averaging origin for reducing the scatter and establishing the tight relationship between massive BHs and bulges (C. Y. Peng 2007; K. Jahnke & A. V. Macciò 2011). This picture is supported by recent findings of a high AGN pair fraction (e.g., M. Perna et al. 2023; J. Li et al. 2024) and the prevalence of companion galaxies around BL AGNs (e.g., J. Matthee et al. 2024) at high redshift.

However, the large scatter may partly be due to our inability to constrain how the mass relation extends to the low-$M_{BH}$ end at a fixed $M_*$, as well as the assumption of a linear $M_{BH}-M_*$ relation with a mass-independent Gaussian scatter, which may not be the case (e.g., M. Habouzit et al. 2022). Given the limited sample size, substantial measurement uncertainties, and significant parameter degeneracies, the recovered intrinsic relation is subject to considerable uncertainties such that both the KH13 and G20 relations fall within its $1\sigma$ uncertainty region. Therefore, it is still possible that the observed
<table><caption>**Table 2**
Model Fitting Results in the Form of log $M_{BH} = \alpha \times \log M_{*} + \beta$ with a Gaussian Intrinsic Scatter $\sigma$ and a BL AGN Fraction $p_{ac}$ (Columns (3)–(6)) Derived Under Different Integration Ranges for $p_{ac}$ and $N_{obs}$ (Columns (1) and (2))</caption><thead><tr><th>p<sub>ac</sub> Range<br>(1)</th><th>N<sub>obs</sub> Range<br>(2)</th><th>α<br>(3)</th><th>β<br>(4)</th><th>σ<br>(5)</th><th>p<sub>ac</sub><br>(6)</th><th>Δ log M<sub>BH</sub> (Default)<br>(7)</th><th>Δ log M<sub>BH</sub> (P23)<br>(8)</th></tr></thead><tbody><tr><td>1%–15%</td><td>15–65</td><td>1.02<sup>+0.56</sup><sub>-0.54</sub></td><td>-3.32<sup>+4.56</sup><sub>-5.16</sub></td><td>0.97<sup>+0.52</sup><sub>-0.37</sub></td><td>0.09<sup>+0.04</sup><sub>-0.05</sub></td><td>...</td><td>-1.26</td></tr><tr><td>10%–15%</td><td>15–65</td><td>1.15<sup>+0.50</sup><sub>-0.60</sub></td><td>-4.72<sup>+4.84</sup><sub>-4.61</sub></td><td>1.04<sup>+0.50</sup><sub>-0.37</sub></td><td>0.13<sup>+0.02</sup><sub>-0.02</sub></td><td>-0.20</td><td>-1.46</td></tr><tr><td>1%–15%</td><td>40–65</td><td>0.91<sup>+0.55</sup><sub>-0.48</sub></td><td>-1.99<sup>+4.00</sup><sub>-5.02</sub></td><td>0.90<sup>+0.51</sup><sub>-0.37</sub></td><td>0.09<sup>+0.04</sup><sub>-0.05</sub></td><td>0.32</td><td>-0.94</td></tr><tr><td>10%–15%</td><td>40–65</td><td>1.07<sup>+0.53</sup><sub>-0.55</sub></td><td>-3.65<sup>+4.39</sup><sub>-4.76</sub></td><td>0.97<sup>+0.52</sup><sub>-0.37</sub></td><td>0.13<sup>+0.02</sup><sub>-0.02</sub></td><td>0.13</td><td>-1.13</td></tr></tbody></table>

**Note.** The fitting results under our default model assumptions are shown in bold. The average offset in BH mass at $7.8 < \log M_{*}/M_{\odot} < 10.6$ relative to the default result and that of F. Pacucci et al. (2023) are shown in columns (7) and (8), respectively.

$M_{BH}/M_{*}$ may not deviate significantly from the intrinsic value. Moreover, adjusting the integration range for $p_{ac}$ and $N_{obs}$ affects the fitting result, as certain mass relations cannot reproduce the assumed AGN number counts. As summarized in Table 2, restricting $p_{ac}$ to a lower range (e.g., 1%–15%) or a higher range (e.g., 10%–15%) shifts the best-fit relation upward or downward, respectively. Adjusting the $N_{obs}$ range in a similar manner ($N_{obs} = 15–65$ versus 40–65) leads to the opposite effect. The fact that the two parameters act in opposite directions means that failing to identify some BL AGNs in the current galaxy sample may not significantly impact the best-fit mass relation, as both $p_{ac}$ and $N_{obs}$ increase in such cases (e.g., $p_{ac} = 10%–15\%$ and $N_{obs} = 40–65$).

# 6. Discussion

## 6.1. Comparison with Previous Work

F. Pacucci et al. (2023) conducted a statistical analysis of the same JDEEP sample and inferred an intrinsic mass relation of $\log M_{BH} = 1.06_{-0.09}^{+0.09} \log M_{*} - 2.43_{-0.83}^{+0.83}$ with a scatter of $0.69_{-0.39}^{+0.69}$ dex. Their best-fit relation, shown as a magenta dashed curve in Figure 5, is ~1.2 dex above our default result at $7.8 < \log M_{*}/M_{\odot} < 10.6$ and ~0.9–1.5 dex above all the fitting results shown in Table 2. Consequently, they concluded that the $M_{BH}/M_{*}$ ratio at $z \sim 5$ is intrinsically higher than that of local AGNs in A. E. Reines & M. Volonteri (2015, hereafter RV15) at a $>3\sigma$ confidence level, with an even greater discrepancy when compared to the G20 relation (Figure 5). This led them to infer that the overall high-redshift BH population is overmassive by ~10–100 times compared to their local counterparts, indicating a lack of undermassive BHs at high redshift.

However, it is important to recognize the limitations of these local samples and acknowledge that their BH and stellar masses are not always consistently measured compared to high-$z$ samples. As detailed in the Appendix, the $M_{BH}-M_{*}$ relation of local AGNs is likely biased low (i.e., missing massive BHs, such as those hosted by massive quiescent galaxies in KH13) due to potential mass-dependent $p_{ac}$. The G20 sample, constructed from local galaxies with mostly dynamical BH mass measurements, is highly inhomogeneous and incomplete, especially at the low-mass end. In fact, overmassive BHs with $M_{BH}/M_{*} \sim 10\%$ at $\log M_{*}/M_{\odot} \sim 9$ have also been reported in the nearby universe (e.g., N. J. Secrest et al. 2017; A. Ferré-Mateu et al. 2021; M. Mezcuá et al. 2023), occupying a similar region in the $M_{BH}-M_{*}$ plane as high-redshift AGNs. However, such a parameter space is not covered by the G20 sample. Thus, it is essential to bear in mind that the local relation itself may be biased and affected by the

complicated selection functions when conducting evolutionary studies (A. Schulze & L. Wisotzki 2011; F. Shankar et al. 2016).

Aware of these caveats and the large uncertainties of our best-fit relation, we find that our best-fit relation aligns with both the RV15 and G20 relations within the $1\sigma$ uncertainty region. There is no evidence that the growth of the entire BH population outpaces their host galaxy at early cosmic times. The observed AGNs do systematically lie above the mean local $M_{BH}-M_{*}$ relation due to selection effects, but they do not represent the entire population of AGNs at high redshift.

The discrepancy between our results and those of F. Pacucci et al. (2023) arises mainly from two differences. First, F. Pacucci et al. (2023) adopted a constant detection limit for BH mass ($\log M_{BH, lim}/M_{\odot} \sim 6.2$ at $z \sim 5$), derived from the $H\alpha$ flux limit (~$10^{-19}$ erg s$^{-1}$ cm$^{-2}$) and a BL width of 1000 km s$^{-1}$, via the average virial estimator. In contrast, our selection function is expressed in terms of AGN luminosity ($\log L_{bol, lim}/\text{erg s}^{-1} \sim 44.1$) and BL width, with the detection of a given BH mass depending on their interplay with the underlying ERDF (Equation (5)). Second, F. Pacucci et al. (2023) allowed the model parameters to vary freely during fitting, whereas we restricted the model parameters to the permitted region.

We display the bivariate $M_{BH}-M_{*}$ and $M_{BH}-\lambda_{Edd}$ distributions from both works in the observed plane (i.e., a mock AGN sample incorporating mass uncertainties and the sample selection function) derived from Equation (6) to demonstrate the differences between the two methodologies. The mock AGN samples derived using the selection function and the best-fit intrinsic $M_{BH}-M_{*}$ relation from this work and F. Pacucci et al. (2023) are denoted as Mock-L and Mock-P, respectively. To ensure consistency with the model assumptions in F. Pacucci et al. (2023), we adopt the M. Song et al. (2016) SMF to simulate their mass relation.

The resulting bivariate distribution, with and without selection functions applied, are represented by the orange and blue contours in Figure 6, respectively. Considering a small grid with $8.0 < \log M_{*}/M_{\odot} < 8.1$ and $6.7 < \log M_{BH}/M_{\odot} < 6.8$, all the mock BHs in that bin are detectable in F. Pacucci et al. (2023) as their masses exceed $\log M_{BH, lim}/M_{\odot} \sim 6.2$. However, only those with $\log L_{bol}/\text{erg s}^{-1} > 44.1$ (see the left inset in Figure 6) and FWHM $> 1000$ km s$^{-1}$ can surpass our detection limit and be included in the Mock-L sample. This different treatment has dramatically altered the distribution of model BHs, thereby affecting the evaluation of the likelihood of the observed sample, even though the average detection limit on $M_{BH}$ is somewhat similar in both works.
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Bivariate distribution of the Mock-L (top) and Mock-P (bottom) samples in the observed (i.e., including mass uncertainties) $M_{BH}-\lambda_{Edd}$ and $M_{BH}-M_*$ planes. Mock samples with (blue contours) and without (orange contours) incorporating selection functions (log $L_{bol} > 44.1$ erg s⁻¹ and FWHM > 1000 km s⁻¹) are shown. The observed JDEEP sample is represented by black points. The distribution of $M_{BH}-\lambda_{Edd}$ for mock AGNs with log $L_{bol} > 44.7$ erg s⁻¹ and FWHM > 1100 km s⁻¹, matching the selection function of the J. Matthee et al. (2024) sample (orange points), is shown in the upper right inset. The inset in the upper left panel presents the luminosity distribution of Mock-L AGNs with $8.0 < \log M_*/M_\odot < 8.1$ and $6.7 < \log M_{BH}/M_\odot < 6.8$. Only BHs with log $L_{bol} > 44.1$ erg s⁻¹ and FWHM > 1000 km s⁻¹ (blue shaded region, corresponding to log $\lambda_{Edd} \gtrsim -0.70$) are included in the final Mock-L sample. However, they all remain in the Mock-P sample since their BH masses exceed the detection limit (log $M_{BH} > 6.2$ $M_\odot$) in F. Pacucci et al. (2023).</figcaption></figure>

Notably, the Mock-L sample (blue contours) successfully replicates the distribution of observed AGNs in both the $M_{BH}-\lambda_{Edd}$ and $M_{BH}-M_*$ planes, demonstrating that only BHs above a limit of $L_{bol} \propto M_{BH} \times \lambda_{Edd}$ can be detected. Our results align well with the predictions from the Cosmic Archaeology Tool (CAT) semi-analytic model (A. Trinca et al. 2023), which suggests that AGNs detected in a CEERS-like survey peak within the mass range of $6 < \log M_{BH}/M_\odot < 8$ and $8 < \log M_*/M_\odot < 10$. The slight deviation between data and model contours in the $M_{BH}-\lambda_{Edd}$ plane arises because these

samples are not strictly flux limit selected; instead, they involve complex target prioritization during shutter allocation. Consequently, the luminosity distribution of mock AGNs does not precisely match the observed samples, with the latter having a slightly higher mean.

In the top right panel of Figure 6, we also display the $M_{BH}-\lambda_{Edd}$ distribution of mock AGNs selected with log $L_{bol}/\text{erg s}^{-1} > 44.7$ and FWHM > 1100 km s⁻¹, which matches the flux-limited BL AGN sample in J. Matthee et al. (2024) selected from slitless spectroscopy. With a shallower
luminosity limit, the $M_{BH}-\lambda_{Edd}$ distribution shifts toward higher BH mass. Notably, our model well reproduces the distribution of the J. Matthee et al. (2024) sample, which has a more uniform sample selection criteria. The predicted observable BL AGN number for the J. Matthee et al. (2024) sample is $N_{obs} \sim 18.1_{-7.4}^{+12.4}$ at $4.2 < z < 5.5$, consistent with the actual observed number of 20, thereby validating our assumed ERDF and fitting result.

Conversely, the Mock-P sample vastly overpredicts the number of “detectable” low-$M_{BH}$, low-$\lambda_{Edd}$ AGNs, which are not observed in real samples. Additionally, although the Mock-P sample generally follows the average $M_{BH}-M_*$ relation of the observed data points, as it is optimized for in the fitting, its highest density region (peaking at $7 < \log M_*/M_{\odot} < 8$; see also Section 4.3 in F. Pacucci et al. 2023) deviates significantly from the observed sample. This discrepancy arises because the statistical model in F. Pacucci et al. (2023) ignores the substantial scatter associated with the virial estimator when deriving the BH mass limit (see Equation (2)). Additionally, it overlooks the fact that a low-mass BH with a high $\lambda_{Edd}$ can still surpass the flux limit, while a massive BH with a low $\lambda_{Edd}$ can be undetectable. This population of low-$M_{BH}$, low-$\lambda_{Edd}$ AGNs, hosted by $\log M_*/M_{\odot} < 8$ galaxies, is responsible for the huge number of “detectable” but missing AGNs with $\log M_{BH, lim} < \log M_{BH} < 6.5$ reported in F. Pacucci et al. (2023). Our simulations and comparison results demonstrate that adopting a constant detection limit for $M_{BH}$ is an inadequate representation of the sample selection function. It is imperative to integrate the ERDF self-consistently into the statistical model to accurately assess the impact of flux limits on the detection limit of BH mass.

## 6.2. Implications for JWST Studies of BH Seedings

Our modeling of selection biases suggests that the observed AGNs may only represent the tip of the iceberg of the underlying BH population. The inferred intrinsic mass relation in the low-mass regime implies the potential existence of a significant population of low-mass BHs with $\log M_{BH}/M_{\odot} \lesssim 5$ that remain undetected or unidentified in current JWST surveys.

The presence of these missing low-mass BHs is supported by the observation that only $\lesssim 10\%$ of galaxies at $4 < z < 7$ in JADES and CEERS host relatively luminous AGNs with $L_{bol} \gtrsim 10^{44} \text{erg s}^{-1}$ (e.g., Y. Harikane et al. 2023; R. Maiolino et al. 2023). Approximating the sample-averaged $\lambda_{Edd}$ ($\sim 0.2$) as the time-averaged $\lambda_{Edd}$ for BH growth, we estimate that a seed BH with $M_{BH} \sim 10^5 M_{\odot}$ would need to accrete at $\lambda_{Edd} \gtrsim 0.2$ with a duty cycle of 100% to grow to $M_{BH} \sim 10^7 M_{\odot}$ by $z \sim 4$. Such a high duty cycle implies that if all galaxies with $8 < \log M_*/M_{\odot} < 10$ at $4 < z < 7$ host overmassive BHs following the mass relation in F. Pacucci et al. (2023), the occupation fraction of AGNs with $\log L_{bol}/\text{erg s}^{-1} \sim 44.4$ ($M_{BH} \sim 10^7 M_{\odot}$, $\lambda_{Edd} \sim 0.2$) should approach $\sim 100\%$. This contrasts with the findings of J. Scholtz et al. (2023), which report a luminous AGN (unobscured + obscured) fraction of $\sim 10\%$. Consequently, BHs in the remaining $\sim 90\%$ of galaxies are likely to have significantly lower masses unless these galaxies universally host heavily obscured luminous AGNs missed by the selection in J. Scholtz et al. (2023), or the bulk of BH growth occurs during a super-Eddington phase at $z > 7$. The AGN fraction inferred from deep multiband MIRI surveys reaching a $5\sigma$

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Comparison of the intrinsic $M_{BH}-M_*$ relation (i.e., true masses without uncertainties) with theoretical predictions at $z \sim 5$. The filled blue stars represent the observed JDEEP AGNs. The filled squares, color coded by the average $L_{bol}$ and the corresponding F277W magnitude in each mass bin, represent mock AGNs simulated using Equation (6) based on our default intrinsic mass relation (black solid line). The filled green triangles (top) and circles (bottom) show the evolutionary tracks (color coded by redshift) of heavy ($10^5 M_{\odot}$) and light ($100 M_{\odot}$) seeds based on the $\tau_{fold} = 225$ Myr model in M. T. Scoggins & Z. Haiman (2024), corresponding to BHs accreting at the mean $\lambda_{Edd}$ of JDEEP AGNs. The $M_{BH}-M_*$ relations of the descendants of heavy and light seeds evolved to $z \sim 5$ in the Eddington-limited CAT model are depicted as orange and green contours ($3\sigma$ level), respectively. The predicted distribution of BHs detectable in the JADES Deepest survey (F277W $< 30.5$ mag and $\log L_{bol} \gtrsim 42.6 \text{ erg s}^{-1}$) under idealized scenarios (see the text for details) is shown as blue contours ($1\sigma$ and $3\sigma$ levels).</figcaption></figure>

depth of $\sim 27$ mag (assuming $\log L_{bol} \sim 44.4 \text{ erg s}^{-1}$ and the type 2 AGN template in J. Lyu et al. 2017) would strongly constrain the occupation fraction of luminous obscured AGNs among $z \sim 4-7$ galaxies, thereby shedding light on the prevalence or absence of undermassive BHs.

Missing a significant population of low-mass BHs has important consequences for understanding BH seeding mechanisms and the early coevolution of SMBHs and their host galaxies. In Figure 7, we display the assembly history of heavy seeds ($10^5 M_{\odot}$) and light seeds ($100 M_{\odot}$) that were seeded at $z=25$ and evolved to a final stellar mass of $\log M_*/M_{\odot} \sim 8.3$ at $z \sim 5$ with an e-folding time of $\tau_{fold} = 225$ Myr, based on the semi-analytic model from M. T. Scoggins & Z. Haiman (2024) (S23). The choice of $\tau_{fold}$ ensures that, on average, BHs are accreting at the mean Eddington ratio (0.2) of JDEEP AGNs. The stellar mass is derived from scaling relations with halo mass, following J. H. Wise et al. (2014) and P. Behroozi et al. (2019), and does not include self-regulated feedback between BHs and galaxies. The comparison of the evolutionary tracks of model seed BHs with our intrinsic $M_{BH}-M_*$ relation shows that while the observed overmassive AGNs may originate from heavy seeds (e.g., F. Pacucci et al. 2023; Á. D. Bogdán et al. 2024; D. Narayanan et al. 2024), lighter seeds are responsible for the majority of low-mass BHs ($\log M_*/M_{\odot} \lesssim 5$) that fall below the average mass relation.

Theoretically, the existence of this low-mass population depends on the effective growth of light- and medium-weight
seeds. B. D. Smith et al. (2018) explored the evolutionary track of light seeds within the framework of the Renaissance simulation and found that the growth of light seeds is impeded by their random motion and distribution in the host galaxy. Moreover, although dense gas clumps are abundant in high-redshift mini-haloes, they are either consumed by star formation and destroyed by stellar feedback, or the light seeds are too small to settle into such clumps (see also M. A. Alvarez et al. 2009). A similar conclusion has been reached in the Eddington-limited CAT model (A. Trinca et al. 2022; R. Schneider et al. 2023), where the inefficient growth of light seeds results in a distinct mass gap between the descendants of light and heavy seeds in the $M_{BH}-M_*$ plane, as shown in Figure 7. However, in models that allow episodic super-Eddington accretion (A. Trinca et al. 2022; R. Schneider et al. 2023), or if ubiquitous cold dense gas flows in the early universe can facilitate continuous but chaotic accretion (thereby reducing BH spin and radiative efficiency) at sufficiently high rates (T. Alexander & P. Natarajan 2014; K. Zubovas & A. King 2021), or if early BH growth is driven by frequent mergers (A. K. Bhowmick et al. 2024), even light seeds can rapidly accumulate significant BH masses.

In the case of medium-weight seeds forming in stellar clusters, their growth is expected to be more efficient compared to light seeds. This is because they can consume gas in the stellar clusters at super-Eddington rates (T. Alexander & P. Natarajan 2014), and their dynamics are stabilized by the dense stellar envelope, which facilitates both their growth and their mergers (P. Biernacki et al. 2017). These efficient growth mechanisms result in a continuous and elevated mass relation that aligns with the observed overmassive AGNs (A. Trinca et al. 2022; A. K. Bhowmick et al. 2024).

Interestingly, the average $M_{BH}-M_*$ relation evolved from heavy seeds in the Eddington-limited CAT model, which represents BHs growing from “typical” evolutionary pathways (i.e., no extreme conditions are required to form significant outliers in $M_{BH}/M_*$), aligns well with our best-fit relation within the observed stellar mass range (Figure 7). However, to explain our mass relation with a larger scatter toward both massive and low-mass BHs at a given stellar mass, more efficient and diverse growth mechanisms are needed. This could include a combination of Eddington-limited and super-Eddington accretion onto different seeds, as well as mergers. Observations directly probing the $M_{BH} < 10^5 M_*$ regime are crucial for validating our intrinsic $M_{BH}-M_*$ relation and distinguishing between theoretical models.

In Figure 7, we present the predicted observable $M_{BH}-M_*$ relation for the JADES Deepest survey (D. J. Eisenstein et al. 2023) under idealized conditions at $z \sim 5$ as blue contours (1$\sigma$ and 3$\sigma$ levels). This is achieved by applying the selection function, which represents the detection limit in AGN bolometric luminosity converted from the $10\sigma$ depth in the F277W band (F277W $\lesssim$ 30.5 mag), to our intrinsic $M_{BH}-M_*$ relation. In this experiment, we assume that all AGNs above the detection limit are free from dust extinction and can be effectively identified, with their BH and stellar masses accurately measured. We acknowledge that this prediction assumes that low-mass BHs follow the same ERDF as their more massive counterparts, which may not hold true. Although the JADES Deepest survey can probe luminosities down to $\log L_{bol}/\text{erg s}^{-1} \sim 42.6$ and unveil BHs with $\log M_{BH}/M_\odot \sim 4$ accreting at the Eddington limit, the rarity of such objects (lying outside the 3$\sigma$ blue contour in

Figure 7) suggests that most detectable BHs will likely have $5 < \log M_{BH}/M_\odot < 6$ and $6 < \log M_*/M_\odot < 8$, accreting at $\log \lambda_{Edd} \sim -0.7$. Consequently, identifying light and medium-weight seeds remains an extremely challenging task even with the deepest JWST surveys. Their presence might only be revealed by detecting gravitational wave emissions from seed BH mergers (e.g., A. Sesana et al. 2009; A. Ricarte & P. Natarajan 2018; A. K. Bhowmick et al. 2024).

It is worth noting that detecting a flat BH distribution at $5 < \log M_{BH}/M_\odot < 6$ and $6 < \log M_*/M_\odot < 8$ in JADES Deepest does not necessarily verify the prediction from the Eddington-limited CAT model at the low-$M_*$ end, which is characterized by a gap and flattening at $\log M_{BH} \sim 5$ (orange and green contours in Figure 7). This is because, due to the finite detection limit, the observable BHs predicted by our intrinsic mass relation—which differs significantly from the CAT model at $\log M_*/M_\odot < 8$—exhibit a nearly identical flat distribution in the $M_{BH}-M_*$ plane (see the 1$\sigma$ blue contour in Figure 7).

In summary, caution is warranted when inferring the dominant formation channel of high-redshift SMBHs from biased samples and comparing observations with theoretical models, as “undermassive” AGNs originating from lighter seeds remain largely undetectable (M. Volonteri et al. 2023).

## 7. Conclusions

In this paper, we investigated the impact of selection biases (i.e., finite detection limits and requirements for detecting BLs) and measurement uncertainties on the observed $M_{BH}-M_*$ relation at high redshift. We demonstrated that these factors can largely account for the observed flattening and enhancement in the ratio of $M_{BH}/M_*$ compared to local values for both JWST-discovered low-luminosity AGNs and the most luminous quasars at $4 < z < 7$ (Section 5.2).

Moreover, by forward modeling the observed low-luminosity AGN sample, we inferred an intrinsic, bias-corrected mass relation in the low-mass regime, parameterized as $\log M_{BH} = 1.02^{+0.56}_{-0.54} \log M_* - 3.32^{+4.56}_{-5.16}$ with a scatter of $0.97^{+0.52}_{-0.37}$. In contrast to F. Pacucci et al. (2023), who concluded that the high-redshift BH population is intrinsically overmassive by $\sim 10-100$ times compared to local counterparts (i.e., lacking undermassive BHs at high-$z$), our careful treatment of selection bias results in an intrinsic $M_{BH}/M_* \sim 1.2$ dex below the observed, biased sample and aligns with local values (Sections 5.3 and 6.1). Our results suggest that the observed AGNs only populate the upper envelope of the underlying BH population in the low-$M_*$ regime, even in deep JWST fields. This highlights the importance of properly considering the scatter of the mass relation when investigating its potential evolution.

The inferred intrinsic mass relation suggests that a significant population of low-mass BHs ($M_{BH} < 10^5 M_\odot$), potentially formed from lighter seeds, might be missing even in the deepest JWST surveys (Section 6.2). Meanwhile, the elevated scatter in the high-redshift $M_{BH}-M_*$ relation in the low-$M_*$ end, compared to the KH13 relation and massive quasars at $z \lesssim 2$ (J. Li et al. 2021a; T. S. Tanaka et al. 2025), could indicate that mergers have played a key role in establishing the close connections between massive BHs and their bulges in the local universe (C. Y. Peng 2007; K. Jahnke & A. V. Macciò 2011).

However, inferring the intrinsic mass relation from the observed, biased sample is a challenging task fraught with
substantial uncertainties. With current data, we cannot rule out the possibility that the observed $M_{BH}/M_\star$ may not deviate significantly from the intrinsic value (Figure 5). Extending the current “effective” luminosity limit to the theoretical minimum achievable in JWST surveys is essential for obtaining more precise constraints on the low-mass BH population. The secondary BHs in the three dual AGN candidates reported by R. Maiolino et al. (2023), with an average $M_{BH} \sim 1.3$ dex smaller than the BL AGN sample analyzed here, support the potential existence of an abundant population of lower-mass BHs accessible with JWST. Significant progress has also been made in identifying type 2 AGN candidates at $z > 4$ with luminosities $\sim 10$ times lower than those of BL AGNs, primarily through the detection of high-ionization lines (J. Scholtz et al. 2023). Assuming that type 2 AGNs with $L_{bol} < 10^{44} \text{erg s}^{-1}$ in J. Scholtz et al. (2023) follow the same accretion rate distribution as type 1 s, their inferred $M_{BH}/M_\star$ is $\sim 0.3\%$, consistent with the local value. Establishing alternative methods for estimating BH masses in these low-mass AGNs—such as narrow line ratios or variability (D. Baron & B. Ménard 2019; M. Sun 2023)—are imperative, as their BLs are obscured or narrowed, making them indistinguishable from those powered by star formation processes.

Besides, our bias correction approach inevitably depends on the assumed forms of the intrinsic mass relation and the underlying distribution functions. Obtaining independent measurements of the ERDF, AGN fraction, and duty cycle (e.g., through clustering measurements) (Y. Shen et al. 2007) in the low-mass regime, as well as improving theoretical predictions on the assembly history of BH seeds, will be crucial for refining our model assumptions, narrowing down the fitting parameter space, and reducing posterior uncertainties.

On the other hand, the JDEEP and SHELLQs surveys have begun to probe BHs positioned on the local relation at the massive stellar mass end. However, massive galaxies are rare in pencil-beam surveys, so the current sample size is insufficient for statistical analysis. Expanding deep JWST surveys to cover a significantly larger area with a uniform sample selection (e.g., Y. Shen et al. 2024b), or extending JWST observations of SHELLQs quasars to even fainter systems would provide the opportunity to directly probe undermassive BHs in massive galaxies at early cosmic times (J. Li et al. 2025). This would allow for a more robust constraint on the intrinsic $M_{BH}-M_\star$ relation at the massive end that mitigates the impact of selection bias.

## Acknowledgments

We thank the referee for very useful comments that helped improve the paper. We thank Raffaella Schneider and Alessandro Trinca for providing us with the model predictions from CAT. J.S. is supported by JSPS KAKENHI (JP22H01262) and the World Premier International Research Center Initiative (WPI), MEXT, Japan. This work was supported by the JSPS Core-to-Core Program (grant No. JPJSCCA20210003). This work is partially supported by JWST programs JWST-GO-02057 and JWST-AR-03038 through a grant from the STScI under NASA contract NAS5-03127. This work was supported in part by NSF grant AST-2006176.

# Appendix

## Impact of a BH Mass-dependent BL AGN Fraction

A. Schulze & L. Wisotzki (2011) first illustrated that a mass-dependent $p_{ac}$ induces a further bias to the observed mass relation when investigating its cosmic evolution. This arises from the fact that high-redshift observations can only detect BHs and measure their masses when they are activated as BL AGNs. The activation efficiency and duration, and hence $p_{ac}$, depends on the complex interplay between the self-regulated fueling and feedback of BHs with their host galaxy and the gas reservoir, and is likely to evolve strongly with redshift as well as BH and stellar masses. Consequently, the mass relation traced by BL AGNs is *not* equivalent to that of the entire BH population. A constant $M_{BH}-M_\star$ relation coupled with a redshift- and mass-dependent active fraction could result in an evolving mass relation for AGNs (A. Schulze & L. Wisotzki 2011; M. Volonteri & D. P. Stark 2011).

In this Appendix, we investigate the impact of a BH mass-dependent $p_{ac}$ on the observed mass relation. Given the lack of observational constraint at high redshift, we simply express the mass dependency as $p_{ac}(m) \propto (0.1m - 0.4)^\gamma$ (arbitrary normalization, which will be canceled out in Equation (7); negative values are set to 0), where $\gamma$ can take values of 1 and 5 (Figure 8(a)). This somewhat arbitrary form is chosen to mimic the increased fueling efficiency and duty cycle of massive BHs at high redshift (R. Aversa et al. 2015; Q. Zhu et al. 2022), which is required to explain their large BH masses assembled within the first billion years of the Universe.

Following the procedure outlined in Section 5.1 and Figure 2, we present in Figure 8(b) the predicted observed (i.e., incorporating measurement uncertainties) mass relation for the entire $z \sim 5$ BL AGN population (i.e., including those below the luminosity limit) by incorporating $\Omega(m) = p_{ac}(m)$ into Equation (6), assuming that the underlying BH population follows the KH13 relation. Both a constant and a BH mass-dependent $p_{ac}$ ($\gamma = 1$ and 5) are displayed for comparison. It can be seen that the predicted mass relation for BL AGNs is biased upwards when a mass-dependent $p_{ac}$ is adopted, even without invoking any selection on AGN luminosity. This is caused by assigning more weight to massive BHs through $p_{ac}$. The effect is particularly pronounced for the $\gamma = 5$ model, where $p_{ac}$ increases dramatically with increasing $M_{BH}$.

We then incorporate the selection function of the JDEEP sample to predict the observed mass relation for detectable AGNs (Figure 8(c)). Interestingly, the resulting mass relation of a constant and a mass-dependent $p_{ac}$ model becomes similar. The median difference of the two mass relations is 0.05 dex for $\gamma = 1$ and 0.25 dex for $\gamma = 5$. This is because the $M_{BH}$ distribution for BHs surpassing the selection criteria is limited to a relatively narrow range (spanning $\sim 1$ dex for the highest density region), where $p_{ac}$ does not vary strongly with the BH mass.

The situation might be reversed in the local universe, as the accretion onto the most massive BHs has almost ceased, such that their BH mass is inferred through dynamic modeling (KH13). In Figure 9, we investigate the impact of a redshift- and BH mass-dependent $p_{ac}$ on the observed mass relation for BL AGNs at $z \sim 5$ and $z \sim 0$, again under the assumption that the underlying BH population intrinsically follows the KH13 relation. We adopt the $\gamma = 5$ model for $z \sim 5$ AGNs. The $p_{ac}(m)$ in the local universe is adopted from T. T. Ananna et al. (2022) based on Swift/BAT AGNs, which can be
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Impact of a BH mass-dependent BL AGN fraction ($p_{ac}$) on the observed (i.e., incorporating mass uncertainties with $\sigma_m = 0.35$ and $\sigma_s = 0.45$) mass relation at $z \sim 5$. Panel (a): $p_{ac}$ as a function of BH mass in the form of $p_{ac} \propto (0.1\,m - 0.4)^\gamma$ with $\gamma = 1$ and $\gamma = 5$. Panel (b): the predicted observed mass relation for the entire BL AGN population (i.e., including those below the flux limits) at $z \sim 5$ based on a constant (orange) and a mass-dependent $p_{ac}$ (blue), assuming that the underlying BH population follows the KH13 relation. The bivariate $M_{BH}-M_*$ distribution is shown as contours, while the mean relation is represented by dashed curves. Panel (c): similar to panel (b) but for detectable BL AGNs with $\log L_{bol} > 44.1$ erg s$^{-1}$ and FWHM > 1000 km s$^{-1}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 9.** The distinct impact of a redshift- and BH mass-dependent BL AGN fraction (panel (a)) on the observed mass relation at $z \sim 5$ and $z \sim 0$. The underlying BH population at both redshifts is assumed to follow the KH13 relation (black dashed line). The blue curves and contours mirror the $\gamma = 5$ model in Figure 8, showcasing the predicted mass relation for BL AGNs at $z \sim 5$ after incorporating the mass-dependent $p_{ac}$ (panel (b)) together with the luminosity and line width cut (panel (c)) to the KH13 relation. The orange curves and contours show the predicted mass relation for $z \sim 0$ BL AGNs by replacing the $\gamma = 5p_{ac}$ model with the Schechter model (orange curve in panel (a)) in T. T. Ananna et al. (2022), and adjusting the selection function to $L_{bol} > 10^{42}\,\text{erg s}^{-1}$ and FWHM > 500 km s$^{-1}$ for local AGNs.</figcaption></figure>

approximated by a Schechter function that decreases strongly at the massive end (Figure 9(a)). As shown in Figure 9(b), the predicted mass relation for local BL AGNs is biased toward

low-mass BHs based on the Schechter $p_{ac}$. Its discrepancy with high-redshift AGNs becomes even more prominent when the selection function, represented by $\log L_{bol}/\text{erg s}^{-1} > 44.1$ and
FWHM > 1000 km s⁻¹ at z ~ 5, and log L<sub>bol</sub>/erg s⁻¹ > 42.0 and FWHM > 500 km s⁻¹ at z ~ 0 (RV15), is incorporated (Figure 9(c)). Therefore, at least part of the offset between high-z and local AGNs seen in F. Pacucci et al. (2023) and M. A. Stone et al. (2024) can be attributed to the mass-dependent p<sub>ac</sub>.

We emphasize that the exact form of p<sub>ac</sub>(m) is still highly uncertain, and thus our results based on two extreme cases of p<sub>ac</sub> at z ~ 0 and z ~ 5 should be treated as illustrative only. Nevertheless, our analysis underscores that BL AGNs are biased tracers of BHs, which must be kept in mind when adopting AGNs as the local baseline to study the evolution of the mass relation.

ORCID iDs

Junyao Li https://orcid.org/0000-0002-1605-915X

John D. Silverman https://orcid.org/0000-0002-0000-6977

Yue Shen https://orcid.org/0000-0003-1659-7035

Marta Volonteri https://orcid.org/0000-0002-3216-1322

Knut Jahnke https://orcid.org/0000-0003-3804-2137

Ming-Yang Zhuang https://orcid.org/0000-0001-5105-2837

Matthew T. Scoggins https://orcid.org/0000-0002-0748-9115

Xuheng Ding https://orcid.org/0000-0001-8917-2148

Yuichi Harikane https://orcid.org/0000-0002-6047-430X

Masafusa Onoue https://orcid.org/0000-0003-2984-6803

Takumi S. Tanaka https://orcid.org/0009-0003-4742-7060

References

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

Aihara, H., AlSayyad, Y., Ando, M., et al. 2022, PASJ, 74, 247

Alexander, T., & Natarajan, P. 2014, Sci, 345, 1330

Alvarez, M. A., Wise, J. H., & Abel, T. 2009, ApJL, 701, L133

Ananna, T. T., Weigel, A. K., Trakhtenbrot, B., et al. 2022, ApJS, 261, 9

Aversa, R., Lapi, A., de Zotti, G., Shankar, F., & Danese, L. 2015, ApJ, 810, 74

Baggen, J. F. W., van Dokkum, P., Brammer, G., et al. 2024, ApJL, 977, L13

Baggen, J. F. W., van Dokkum, P., Labbé, I., et al. 2023, ApJL, 955, L12

Baron, D., & Ménard, B. 2019, MNRAS, 487, 3404

Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, MNRAS, 370, 289

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

Bhowmick, A. K., Blecha, L., Torrey, P., et al. 2024, arXiv:2402.03626

Biernacki, P., Teysnier, R., & Bleuler, A. 2017, MNRAS, 469, 295

Bogdán, Á., Goulding, A. D., Natarajan, P., et al. 2024, NatAs, 8, 126

Chabrier, G. 2003, PASP, 115, 763

Ding, X., Onoue, M., Silverman, J. D., et al. 2023, Natur, 621, 51

Ding, X., Silverman, J. D., Treu, T., et al. 2022, ApJ, 933, 132

Du, P., Lu, K.-X., Zhang, Z.-X., et al. 2016, ApJ, 825, 126

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

Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv:2306.02465

Fan, X., Bañados, E., & Simcoe, R. A. 2023, ARA&A, 61, 373

Ferré-Mateu, A., Mezcua, M., & Barrows, R. S. 2021, MNRAS, 506, 4702

Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJL, 946, L13

Fryer, C. L., Woosley, S. E., & Heger, A. 2001, ApJ, 550, 372

Giménez-Arteaga, C., Fujimoto, S., Valentino, F., et al. 2024, A&A, 686, A63

Giménez-Arteaga, C., Oesch, P. A., Brammer, G. B., et al. 2023, ApJ, 948, 126

Goulding, A. D., Greene, J. E., Setton, D. J., et al. 2023, ApJL, 955, L24

Greene, J. E., & Ho, L. C. 2005, ApJ, 630, 122

Greene, J. E., Labbe, I., Goulding, A. D., et al. 2024, ApJ, 964, 39

Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257

Habouzit, M., Onoue, M., Bañados, E., et al. 2022, MNRAS, 511, 3751

Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, ApJ, 959, 39

Hirschmann, M., Khochfar, S., Burkert, A., et al. 2010, MNRAS, 407, 1016

Jahnke, K., & Macciò, A. V. 2011, ApJ, 734, 92

Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, ApJL, 954, L4

Kokorev, V., Chisholm, J., Endsley, R., et al. 2024, ApJ, 975, 178

Kokorev, V., Fujimoto, S., Labbe, I., et al. 2023, ApJL, 957, L7

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

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

Li, J., Shen, Y., & Zhuang, M.-Y. 2025, arXiv:2502.05048

Li, J., Silverman, J. D., Ding, X., et al. 2021a, ApJ, 922, 142

Li, J., Silverman, J. D., Ding, X., et al. 2021b, ApJ, 918, 22

Li, J., Silverman, J. D., Izumi, T., et al. 2022, ApJL, 931, L11

Li, J., Zhuang, M.-Y., Shen, Y., et al. 2024, arXiv:2405.14980

Li, J. I. H., Shen, Y., Ho, L. C., et al. 2023, ApJ, 954, 173

Lyu, J., Rieke, G. H., & Shi, Y. 2017, ApJ, 835, 257

Madau, P., & Rees, M. J. 2001, ApJL, 551, L27

Maiolino, R., Scholtz, T., Curtis-Lake, E., et al. 2023, arXiv:2308.01230

Matsuoka, Y., Onoue, M., Kashikawa, N., et al. 2016, ApJ, 828, 26

Matthee, J., Naidu, R. P., Brammer, G., et al. 2024, ApJ, 963, 129

Mezcua, M., Siudek, M., Suh, H., et al. 2023, ApJL, 943, L5

Narayanan, D., Lower, S., Torrey, P., et al. 2024, ApJ, 961, 73

Natarajan, P., Pacucci, F., Ferrara, A., et al. 2017, ApJ, 838, 117

Navarro-Carrera, R., Rinaldi, P., Caputi, K. I., et al. 2024, ApJ, 961, 207

Onoue, M., Kashikawa, N., Matsuoka, Y., et al. 2019, ApJ, 880, 77

Pacucci, F., Nguyen, B., Carniani, S., Maiolino, R., & Fan, X. 2023, ApJL, 957, L3

Peng, C. Y. 2007, ApJ, 671, 1098

Perna, M., Arribas, S., Lamperti, I., et al. 2023, arXiv:2310.03067

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

Ricarte, A., & Natarajan, P. 2018, MNRAS, 481, 3278

Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al. 2006, ApJS, 166, 470

Schneider, R., Valiante, R., Trinca, A., et al. 2018, MNRAS, 526, 3250

Scholtz, J., Maiolino, R., D'Eugenio, F., et al. 2023, arXiv:2311.18731

Schulze, A., & Wisotzki, L. 2011, A&A, 535, A87

Schulze, A., & Wisotzki, L. 2014, MNRAS, 438, 3422

Scoggins, M. T., & Haiman, Z. 2024, MNRAS, 531, 4584

Scoggins, M. T., Haiman, Z., & Wise, J. H. 2023, MNRAS, 519, 2155

Secrest, N. J., Schmitt, H. R., Blecha, L., Rothberg, B., & Fischer, J. 2017, ApJ, 836, 183

Sesana, A., Gair, J., Mandel, I., & Vecchio, A. 2009, ApJL, 698, L129

Shankar, F., Bernardi, M., Sheth, R. K., et al. 2016, MNRAS, 460, 3119

Shen, Y. 2013, BASI, 41, 61

Shen, Y., Greene, J. E., Ho, L. C., et al. 2015, ApJ, 805, 96

Shen, Y., Grier, C. J., Horne, K., et al. 2024a, ApJS, 272, 26

Shen, Y., Strauss, M. A., Oguri, M., et al. 2007, AJ, 133, 2222

Shen, Y., Zhuang, M.-Y., Li, J., et al. 2024b, arXiv:2408.12713

Smith, B. D., Regan, J. A., Downes, T. P., et al. 2018, MNRAS, 480, 3762

Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5

Steinhardt, C. L., Sneppen, A., Hensley, H., et al. 2022, ApJ, 934, 22

Stern, J., & Laor, A. 2012, MNRAS, 426, 2703

Stone, M. A., Lyu, J., Rieke, G. H., Alberts, S., & Hainline, K. N. 2024, ApJ, 964, 90

Sun, M. 2023, MNRAS, 521, 2954

Sun, M., Trump, J. R., Brandt, W. N., et al. 2015, ApJ, 802, 14

Tacchella, S., Johnson, B. D., Robertson, B. E., et al. 2023, MNRAS, 522, 6236

Tanaka, T. S., Silverman, J. D., Ding, X., et al. 2025, ApJ, 979, 215

Topping, M. W., Stark, D. P., Endsley, R., et al. 2022, MNRAS, 516, 975

Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110

Treu, T., Woo, J.-H., Malkan, M. A., & Blandford, R. D. 2007, ApJ, 667, 117

Trinca, A., Schneider, R., Maiolino, R., et al. 2023, MNRAS, 519, 4753

Trinca, A., Schneider, R., Valiante, R., et al. 2022, MNRAS, 511, 616

Übler, H., Maiolino, R., Pérez-González, P. G., et al. 2024, MNRAS, 531, 355

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

Visbal, E., & Haiman, Z. 2018, ApJL, 865, L9

Vito, F., Brandt, W. N., Yang, G., et al. 2018, MNRAS, 473, 2378

Volonteri, M. 2010, A&ARv, 18, 279

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

Volonteri, M., & Natarajan, P. 2009, MNRAS, 400, 1911

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

Whitler, L., Stark, D. P., Endsley, R., et al. 2023, MNRAS, 519, 5859

Wise, J. H., Demchenko, V. G., Halicek, M. T., et al. 2014, MNRAS, 442, 2560

Wu, J., Shen, Y., Jiang, L., et al. 2022, MNRAS, 517, 2659

Yue, M., Eilers, A.-C., Simcoe, R. A., et al. 2024, ApJ, 966, 176

Zhang, H., Behroozi, P., Volonteri, M., et al. 2023, MNRAS, 523, L69

Zhu, Q., Li, Y., Li, Y., et al. 2022, MNRAS, 514, 5583

Zhuang, M.-Y., & Ho, L. C. 2023, NatAs, 7, 1376

Zhuang, M.-Y., & Shen, Y. 2024, ApJ, 962, 139

Zubovas, K., & King, A. 2021, MNRAS, 501, 4289