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

# Glimmers in the Cosmic Dawn. II. A Variability Census of Supermassive Black Holes across the Universe*

Vieri Cammelli<sup>1,2,3,4,5</sup>,<sup>6</sup> , Jonathan C. Tan<sup>4,6</sup>,<sup>7</sup> , Alice R. Young<sup>7</sup>,<sup>8</sup> , Matthew J. Hayes<sup>7</sup>,<sup>9</sup> , Jasbir Singh<sup>8</sup>,<sup>10</sup> , Richard S. Ellis<sup>9</sup>,<sup>11</sup> , Aayush Saxena<sup>10</sup>,<sup>12</sup> , Nicolas Laporte<sup>11</sup>,<sup>13</sup> , Pierluigi Monaco<sup>2,3,5,12</sup>,<sup>15</sup> , and Benjamin W. Keller<sup>13</sup> <sup>14</sup>

<sup>1</sup> Department of Physics, Informatics & Mathematics, University of Modena & Reggio Emilia, via G. Campi 213/A, 41125, Modena, Italy;
vieri.cammelli@unimore.it

<sup>2</sup> Dipartimento di Fisica, Sezione di Astronomia, Università degli Studi di Trieste, via G.B. Tiepolo 11, I-34131, Trieste, Italy

<sup>3</sup> INAF-Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, I-34131, Trieste, Italy

<sup>4</sup> Department of Space, Earth & Environment, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden

<sup>5</sup> IFPU—Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34151 Trieste, Italy

<sup>6</sup> Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA

<sup>7</sup> Stockholm University, Department of Astronomy and Oskar Klein Centre for Cosmoparticle Physics, AlbaNova University Centre, SE-10691 Stockholm, Sweden
<sup>8</sup> INAF—Astronomical Observatory of Brera, via Brera 28, I-20121 Milan, Italy

<sup>9</sup> Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

<sup>10</sup> Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

<sup>11</sup> Aix Marseille Université, CNRS, CNES, LAM (Laboratoire d'Astrophysique de Marseille), UMR 7326, 13388 Marseille, France

<sup>12</sup> INFN, Sezione di Trieste, Via Valerio 2, I-34127 Trieste, Italy

<sup>13</sup> Department of Physics and Materials Science, University of Memphis, 3720 Alumni Avenue, Memphis, TN 38152, USA

Received 2025 February 1; revised 2025 July 17; accepted 2025 July 29; published 2025 September 23

## Abstract

Understanding the origin and evolution of supermassive black holes (SMBHs) stands as one of the most important challenges in astrophysics and cosmology, with little current theoretical consensus. Improved observational constraints on the cosmological evolution of SMBH demographics are needed. Here we report results of a search via photometric variability for SMBHs appearing as active galactic nuclei (AGN) in the cosmological volume defined by the Hubble Ultra Deep Field. This work includes particular focus on a new observation carried out in 2023 with the Hubble Space Telescope using the WFC3/IR/F140W, which is compared directly to equivalent data taken 11 yr earlier in 2012. Two earlier pairs of observations from 2009 to 2012 with WFC3/IR/F105W and WFC3/IR/F160W are also analyzed. We identify 521, 188, and 109 AGN candidates as nuclear sources that exhibit photometric variability at a level of 2σ, 2.5σ, and 3σ, respectively, in at least one filter. This sample includes 13, 3, and 2 AGN candidates at redshifts $z > 6$, when the Universe was $\lesssim 900$ Myr old. After variability and luminosity function (down to $M_{UV} = -17$ mag) completeness corrections, we estimate the comoving number density of SMBHs, $n_{SMBH}(z)$. At $z \gtrsim 6$, $n_{SMBH} \gtrsim 6 \times 10^{-3}$ cMpc$^{-3}$. At low $z$ our observations are sensitive to AGN fainter than $M_{UV} = -17$ mag, and we estimate $n_{SMBH} \gtrsim 10^{-2}$ cMpc$^{-3}$. We discuss how these results place strong constraints on a variety of SMBH seeding theories.

*Unified Astronomy Thesaurus concepts: High-redshift galaxies (734); Galaxies (573); Active galactic nuclei (16); Supermassive black holes (1663)*

## 1. Introduction

Understanding the origin of supermassive black holes (SMBHs) is one of the most important unsolved problems of astrophysics and cosmology. These SMBHs reside in the nuclei of most large galaxies and are typically detected when undergoing accretion and appearing as active galactic nuclei (AGN). Feedback from SMBHs/AGN may play crucial roles in galaxy assembly (e.g., J. Kormendy & L. C. Ho 2013; M. den Brok et al. 2015; A. W. Graham 2016).

The overall abundance, e.g., as measured by their comoving number density $n_{SMBH}$, of SMBHs across cosmic time is of

crucial importance for understanding their formation and evolution. An approximate estimate of SMBH abundance in the local Universe is $n_{SMBH} \sim 4.6 \times 10^{-3}$ cMpc$^{-3}$ (N. Banik et al. 2019), which assumes that every $\sim L_{*}/3$ galaxy hosts an SMBH, where $L_{*}$ is the characteristic luminosity. Integrating the SMBH mass function of M. Vika et al. (2009) yields a similar estimate of $\sim 8.8 \times 10^{-3}$ cMpc$^{-3}$. However, an undetected population of fainter SMBHs would increase these numbers. The lowest-mass examples of known SMBHs have masses just below $\sim 10^{5} M_{\odot}$. For example, the study of M. Häberle et al. (2024) implies a black hole mass in the range $\sim (2-5) \times 10^{4} M_{\odot}$ in $\omega$ Cen, which is the stripped nucleus of a dwarf galaxy accreted by the Milky Way. The mass function of SMBHs derived from the modeling of tidal disruption events by A. Mummery & S. van Velzen (2025) extends down to $\sim 4 \times 10^{4} M_{\odot}$. Furthermore, from this work there are hints of a turnover in the mass function on scales of $\sim 4 \times 10^{5} M_{\odot}$. This, together with the general lack of good examples of intermediate-mass black holes (IMBHs) in the mass range of $\sim 10^{3}-10^{4} M_{\odot}$, suggests that the initial seed masses for SMBHs may be already in or near the “supermassive” regime of $\sim 10^{5} M_{\odot}$.

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

<footnote>* This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with programs 1563, 12498, and 17073.</footnote>
SMBHs with masses ~$10^9 M_\odot$ have been found in the early Universe at $z \ge 7$ (e.g., A. Mortlock et al. 2013; E. Bañados et al. 2018; F. Wang et al. 2021). The James Webb Space Telescope (JWST) has accelerated claims for high-$z$ SMBHs, including lower-mass (~$10^7–10^8 M_\odot$) examples out to higher redshifts ($z \sim 10$; e.g., J. E. Greene et al. 2024; Y. Harikane et al. 2023b; J. Matthee et al. 2024; R. Maiolino et al. 2024; A. Bogdan et al. 2024; V. Kokorev et al. 2024). Assuming Eddington-limited growth, these results also imply that at least some SMBHs form early with seed masses in the supermassive regime with ~$10^5 M_\odot$ (e.g., F. Wang et al. 2021).

There are several proposed theories for how SMBHs form. Much attention has focused on monolithic “direct collapse” (DC) in metal-free, irradiated, relatively massive (~$10^8 M_\odot$), atomically cooled halos (e.g., M. G. Haehnelt & M. J. Rees 1993; M. C. Begelman et al. 2006; S. Chon et al. 2016; J. H. Wise et al. 2019). However, the conditions required for this mechanism appear to be quite rare. For example, S. Chon et al. (2016) found only two candidate DC SMBHs in their simulation of a ~$(30 \text{ cMpc})^3$ volume, implying $n_{\text{SMBH}} \sim 9 \times 10^{-5} \text{ cMpc}^{-3}$. In the radiation-hydro simulation of J. H. Wise et al. (2019), an even smaller global number density of DC SMBH seeds was inferred of $n_{\text{SMBH}} \sim 10^{-7}$ to $10^{-6} \text{ cMpc}^{-3}$. A variant of the DC scenario in which strong turbulence induced by fast converging flows supports gas in the halo against fragmentation has been proposed by M. A. Latif et al. (2022). However, again, the conditions for such halos appear to be very rare, yielding $n_{\text{SMBH}} \lesssim 8 \times 10^{-7} \text{ cMpc}^{-3}$.

An alternative theory of SMBH formation is based on runaway mergers of stars in dense star clusters (e.g., M. Freitag et al. 2006; D. R. G. Schleicher et al. 2023). However, the stellar densities required to produce stars that would be progenitors of even IMBHs, i.e., with ~$10^3 M_\odot$, are almost never seen in local examples of young, massive star clusters (J. C. Tan et al. 2014). Thus, while it is difficult to predict the occurrence rate of this mechanism in cosmological volumes, it appears to require very rare, specialized conditions. Furthermore, there is very little evidence for clear examples of IMBHs with masses ~$10^3–10^4 M_\odot$, which would be predicted to be much more common than the SMBHs that form via this mechanism.

Given these theoretical uncertainties, many implementations of SMBH seeding in numerical simulations have utilized simple threshold conditions. For example, in the Illustris project simulations (M. Vogelsberger et al. 2014) a halo mass threshold (HMT) condition has been used (based on, e.g., D. Sijacki et al. 2007; T. D. Matteo et al. 2008): when a dark matter halo exceeds $7 \times 10^{10} M_\odot$, it is seeded with an SMBH of mass ~$10^5–10^6 M_\odot$. A natural feature of these models is that SMBH formation occurs relatively late in the Universe, since it takes time for these massive halos to develop. For example, in the fiducial Illustris model, the comoving number density increases from $n_{\text{SMBH}} \sim 10^{-5} \text{ cMpc}^{-3}$ at $z = 10$ to ~$4 \times 10^{-3} \text{ cMpc}^{-3}$ by $z = 5$. It continues to rise toward lower redshifts, asymptotically approaching ~$2 \times 10^{-2} \text{ cMpc}^{-3}$ by $z = 0$.

In another approach, SMBH seeding recipes based on threshold conditions in the gas have been implemented. For instance, in the Romulus25 simulation suite (M. Tremmel et al. 2017) of $(25 \text{ cMpc})^3$ volumes SMBH seeds of $10^6 M_\odot$ are created in gas that has relatively low metallicity ($Z < 3 \times 10^{-4}$), a H nuclei number density $n_H > 3 \text{ cm}^{-3}$, and a temperature between 9500 and 10,000 K. These

conditions are designed to seed SMBHs in gas that is collapsing relatively quickly, but still having low cooling rates. However, it should be noted that these thresholds, especially of density, describe conditions that are very far from those needed to resolve the detailed processes expected to be occurring in supermassive star and/or dense star cluster formation leading to SMBHs.

An alternative theoretical model of SMBH seeding in cosmological volumes is based on the formation of supermassive Population III.1 stars in locally isolated dark matter minihalos (N. Banik et al. 2019; J. Singh et al. 2023; V. Cammelli et al. 2025; see J. C. Tan et al. 2024, for a review). It relies on the physical mechanism of dark matter annihilation heating to change the structure of the protostar (D. Spolyar et al. 2008; J. C. Tan 2008; K. Freese et al. 2009, 2010; A. Natarajan et al. 2009; T. Rindler-Daller et al. 2015). In particular, if the protostar can be kept in a large, swollen state relative to that of the zero-age main-sequence structure, then it may avoid the strong ionizing feedback that acts to limit the growth of “standard” Population III stars (J. C. Tan & C. F. McKee 2004; C. F. McKee & J. C. Tan 2008; J. C. Tan et al. 2010; T. Hosokawa et al. 2011; S. Hirano et al. 2014; H. Susa et al. 2014). The Population III.1 model for supermassive star formation requires adiabatic contraction of the dark matter density via monolithic, relatively slow contraction of the baryons in the center of a minihalo, which has been seen in many numerical simulations of the process (T. Abel et al. 2002; V. Bromm et al. 2002; N. Yoshida et al. 2003). The baryons cool via emission of H$_2$ and HD rotational-vibrational transitions, with trace amounts of these species formed in the gas phase, catalyzed by the presence of residual free electrons following cosmic recombination. On the other hand, UV-irradiated metal-free minihalos, i.e., Population III.2 sources, are thought to have elevated production of H$_2$ and HD leading to fragmentation to small clusters of lower-mass, ~$10 M_\odot$ stars (T. H. Greif & V. Bromm 2006; J. L. Johnson & V. Bromm 2006), which would then not lead to significant adiabatic contraction of the dark matter density and thus a much smaller, likely negligible, impact of WIMP annihilation on the protostellar evolution.

The above distinction between “undisturbed” Population III.1 and irradiated Population III.2 minihalos motivates the main parameter of the Population III.1 model for SMBH seeding, i.e., the isolation distance from previous generations of stars, $d_{\text{iso}}$, that is needed for a minihalo to be a Population III.1 source. N. Banik et al. (2019) modeled Population III.1 sources in a ~$(60 \text{ cMpc})^3$ volume and considered a range of values of $d_{\text{iso}}$ from 10 to 300 kpc (proper distance). They found that a value of $d_{\text{iso}} = 100 \text{ kpc}$ led to a comoving number density of SMBHs of $n_{\text{SMBH}} \simeq 10^{-2} \text{ cMpc}^{-3}$, i.e., consistent with the observational constraints discussed above. A value of $d_{\text{iso}} = 50 \text{ kpc}$ led to an approximately 10 times larger number density of $n_{\text{SMBH}} \simeq 10^{-1} \text{ cMpc}^{-3}$. In both cases, the bulk of the SMBH population formed before $z = 20$ and then maintained a near-constant number density toward lower redshifts. While the N. Banik et al. (2019) study only followed the evolution down to $z = 10$, J. Singh et al. (2023) were able to run the simulation down to $z = 0$, as well as tracking mergers between SMBH-seeded halos. As a consequence of the initially widely separated SMBH locations, i.e., typical comoving separation of 3 Mpc in the $d_{\text{iso}} = 100 \text{ kpc}$ case, mergers only became significant at low redshift, i.e., $z \lesssim 2$. However, even in the $d_{\text{iso}} = 50 \text{ kpc}$ case, $n_{\text{SMBH}}$ only dropped
by about 30% owing to mergers by $z = 0$. Thus, the key prediction of the Population III.1 model of SMBH seeding is a near-constant value of $n_{SMBH}$ at all redshifts up to $z \sim 20$.

This prediction for $n_{SMBH}$ motivates the need for better observational constraints at high redshifts. In this series of papers we examine AGN activity in these early, high-redshift systems ($z > 6$) through photometric variability. Specifically we address this by reimaging the historic deep Hubble Space Telescope (HST) field with HST.

The Hubble Ultra Deep Field (HUDF) is the deepest field for which there is a long history of exquisite HST observations. In order to probe variability to the highest redshifts, near-infrared (NIR) imaging with the Wide Field Camera 3 (WFC3) is essential, and the two most relevant historic data sets are the HUDF09 (GO 11563; PI: Illingworth; 192 orbits) and HUDF12 (GO 12498; PI: Ellis; 30 orbits) imaging campaigns.

In Paper I of this series (M. J. Hayes et al. 2024) we reported first results from a new observation of the HUDF that was made in 2023 in the F140W filter (GO 17073; PI: Hayes; hereafter HUDF23). The images were taken to an equivalent depth to the HUDF12 observation and thus permit a detailed study of variability between the two epochs. An analysis of variability between the HUDF09 and HUDF12 epochs was also carried out. The main result of Paper I was the report of three high-$z$ (i.e., between $z = 6$ and 7) AGN candidates identified via variability and implication of these results for $n_{SMBH}$ and thus SMBH seeding models.

Here in Paper II we present the full variability analysis of the HUDF09, HUDF12, and HUDF23 data sets, including detection of AGN candidates at various thresholds of significance and over the full redshift range of the source population. The plan of the paper is as follows: In Section 2 we introduce the three key data sets used in our study. For each of these we describe the processing pipeline and photometric catalogs. We also discuss our techniques for identifying variable candidates. Section 3 presents our results, including the implications for SMBH seeding models. We present our conclusions in Section 4.

## 2. Data Sets and Methods

Here we review the observation data sets and the methodology used in our study (see also Paper I).

### 2.1. Observations and Data Reduction

The HUDF field was initially observed in the optical with the Advanced Camera for Surveys (ACS; S. V. W. Beckwith et al. 2006), followed by key NIR mode observations with WFC3/IR during 2009–2010 under HUDF09 (R. J. Bouwens et al. 2011), capturing images in three filters (F105W, F125W, and F160W) over 192 orbits. In 2012, the field was reimaged under HUDF12 (R. S. Ellis et al. 2013), significantly deepening the F105W and F160W exposures and adding a fourth filter, F140W, to search for Lyman break galaxies at $z \sim 8$.

To search for photometric variability across all sources within the HUDF IR footprint, we reimaged the field in 2023 September using the F140W filter, replicating the center, field orientation, and depth (30 orbits) of the HUDF12 observation. We processed the F140W image with the calwfc3 pipeline and astrodrizzle (STSCI Development Team 2012)

<table><caption>**Table 1**
Observing Epochs and Times</caption><thead><tr><th>Year</th><th>Filters</th><th>Orbits</th><th>GO# / PI</th></tr></thead><tbody><tr><td rowspan="2">2009–2010</td><td>F105W</td><td>24</td><td>11563/Illingworth</td></tr><tr><td>F160W</td><td>53</td><td>...</td></tr><tr><td rowspan="3">2012</td><td>F105W</td><td>72</td><td>12498/Ellis</td></tr><tr><td>F140W</td><td>30</td><td>...</td></tr><tr><td>F160W</td><td>26</td><td>...</td></tr><tr><td>2023</td><td>F140W</td><td>30</td><td>17073/Hayes</td></tr></tbody></table>

**Note.** This table lists the observing epochs, filters used, orbits, and the corresponding GO numbers and PIs.

software, using High-Level Science Products (HLSP) from the Mikulski Archive for Space Telescopes (MAST) as reference images. Additionally, we reprocessed the F140W image from the UDF12 campaign to confirm that our methods align with the depth of the HLSP image. Concurrently, we independently reprocessed the F105W and F160W images from the 2009 and 2012 epochs to search for variable sources over the shorter, earlier time baseline. This setup enables variability searches across three epochs: the period from 2009 to 2012 is covered by the F105W and F160W filters, while the span from 2012 to 2023 is sampled by the F140W filter alone (see Table 1 for details).

### 2.2. Photometry

For a given filter and time baseline, we run SEtractor (E. Bertin & S. Arnouts 1996) on both epochs. We use an rms map based on the weight map produced by astrodrizzle for each filter, which is itself an inverse variance image based on the input exposures that contributed to each pixel. For each detection, SEtractor outputs the geometric properties of the detection ellipse, specifically WCS and pixel coordinates, and photometric properties of the detected object, namely the flux, the magnitude, and their respective errors estimated accordingly.

Rather than using SEtractor in the standard single mode, we opt for the dual-mode approach. In contrast to Paper I, where we implemented the dual-mode approach using the single epochs themselves as detection images in both directions, here we employ a common detection image for any given filter and epoch. Specifically, we use a combined stacked frame from multiple epochs for the unconvolved F105W, F140W, and F160W images, along with the original HLSP in F125W, all matched to the same WCS. This extremely deep detection image, reaching nearly 300 orbits, provides a robust estimate of galaxy morphologies and the highest-precision centroiding on galaxy nuclei where variability is expected. In addition, this method provides a “global” list of source coordinates and apertures that remains consistent across different runs of SEtractor, ensuring one-to-one correspondence for all extracted photometric catalogs.

Instead of conducting “extended” galaxy photometry, such as using Kron-like or moment-centered apertures as done by SEtractor, we aim to obtain photometry centered on the brightest, unresolved sources within each galaxy determined from the ultradeep stack. To achieve this, we rely on the coordinates of the barycenter (centroid) pixels reported by SEtractor and the background images it generates during
its background subtraction process. At these coordinates, we perform aperture photometry within 4-pixel-diameter (0.26) apertures, applying local background subtraction to exclude local (noncompact) galaxy light.

We then assemble a photometric catalog for each individual filter at each epoch by running SExtractor in the dual mode. As stated above, the dual mode assures us that a possible photometric variability, if any, comes from the same region of the sky, once every frame has been co-aligned with respect to the HLSP. Hence, for each image, we obtain photometry of each galaxy nucleus, where AGN variability is expected. We consider that any off-center variability would likely be due to supernovae (SNe). We defer analysis of such variability to a future paper in this series.

We then compare these local aperture magnitudes in each image relative to the photometric uncertainty of each source. Additionally, the drizzling process causes artificially underestimated uncertainties due to the correlation of signal between adjacent pixels. We account for this effect by multiplying the uncertainties by a correction factor, described below, which is assumed to be constant over the processed area (e.g., S. Casertano et al. 2000; A. S. Fruchter & R. N. Hook 2002). It is worth emphasizing that our photometric measurements do not aim for an accurate estimation in absolute terms, but rather for a comparison between two different epochs in terms of relative variation.

## 2.3. Identifying Variable Sources

Following Paper I, we identify variables in each matched pair of filters using two techniques: comparing the nuclear/ central photometry of galaxies at different epochs and detecting residual sources in pair-subtracted images in any given filter. Regions near the image edges are excluded owing to the dithering pattern causing excess noise, and we focus solely on a central region of the HUDF covering 123" × 139".

### 2.3.1. Photometric Variables

The photometric comparison method follows a similar approach to that of R. O'Brien et al. (2024). We first correct for a systematic offset of approximately 0.01 mag that is observed at all magnitudes, likely due to slightly imperfect zero-points in images taken many years apart, as well as telescope expansion and/or breathing. Under the assumption that the majority of sources in the field will not vary, we calculate the standard deviation of the Δm distribution in 0.5 mag bins, where Δm = m₁ - m₂ is the magnitude difference. The magnitude labels always refer to the time ordering of the different epochs, with “1” being the first visit and “2” being the second. We then propagate the uncertainty on the magnitude difference derived from our photometric catalogs in the following way:

$$ \delta_{\Delta m} = \left( \frac{2.5}{\ln 10} \right) \sqrt{\left(\frac{\delta_{F_1}}{F_1}\right)^2 + \left(\frac{\delta_{F_2}}{F_2}\right)^2}, \quad (1) $$

where δF and F are the flux error and the flux as estimated by SExtractor, respectively.

To calibrate our estimates of the uncertainties, we compare the observed standard deviations of the Δm distributions in magnitude bins of 0.5 mag width over the range from 24 to 31

mag with the estimated mean photometric uncertainties in each bin from Equation (1).

We then compute the total average scale factor that globally adjusts the uncertainty reported in Equation (1) to match the observed Δm standard deviation. For the F105W, F140W, and F160W filters the resulting calibration scale factors are 0.80, 0.84, and 0.79, respectively. We refer to the calibrated uncertainties as δΔm,cal.

It is by taking the ratio of the Δm with the uncertainty calculated in Equation (1) that a given source's significance of variability is ultimately assessed in a two-epoch observation in a given filter. We verify that this empirical approach follows the trend of the calibrated uncertainties by carrying out a linear fit to the logarithm of the calibrated 1σ uncertainty estimates in each 0.5 mag bin (see Figure 1). We consider three thresholds of significance for variability: 2σ, 2.5σ, and 3σ. For a source to be classified as a candidate variable, it must have a σΔm ≡ Δm/δΔm,cal ≥ 2 in at least one filter.

### 2.3.2. Difference Imaging Variables

As a complement to variability identified via aperture photometry, we also generate pair-subtracted images for each filter. From these images we may identify additional variable sources, including sources offset from galactic nuclei, that were not detected as variables via aperture photometry. In addition, for extended sources we are able to assess the morphology of photometric variables. In particular, single AGN are expected to appear as compact, unresolved variable sources. Then, images taken in the same filters were directly subtracted, resulting in six “difference images.” We then applied SExtractor to these difference images, and effective gains were recalculated using the difference of the squared exposure times multiplied by the CCD gain to account for the increased sky noise in the subtracted images, which mimics a shallower integration time. Difference-image-selected sources, as we describe below, are relatively few in number but are considered to have high significance and, for the purposes of counting statistics, are included in the 3σ significance photometric sample.

Finally, to prevent possible image artifacts and other forms of interlopers contaminating our sample, we visually inspect every source identified as variable in both methods. Possible interlopers include moving objects, noisy pixels, saturated sources, and telescope diffraction spikes. SNe that occur between the two epochs (either between 2009 and 2012 or between 2012 and 2023) are removed from the sample since there will be no detection in the first epoch and a bright, obvious detection in the second. For this visual inspection, along with pair-subtracted images, we check the single-epoch frames in all three filters. Furthermore, we rule out the possibility of low-redshift interlopers or artifacts by examining cutouts in the optical and in the infrared by leveraging available observations with HST and JWST, respectively (see Section 3.2 for the used criteria). Sources identified as image artifacts or interlopers are excluded from the variable source sample.

### 2.3.3. Redshift Determination

We cross-correlate each analyzed source with positions from known redshift catalogs, including spectroscopic data from Very Large Telescope/MUSE (R. Bacon et al. 2023), the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top panel: photometric variability of HUDF sources in the F140W filter from 2012 to 2023, Δ*m*<sub>F140W</sub>, vs. mean magnitude, *m*̄<sub>F140W</sub> (gray points). Black circles show the intrinsic 1σ scatter of the data measured in bins of 0.5 mag width. Yellow circles depict the calibrated (by a factor of 0.84) photometric 1σ uncertainties in each bin (see text). The inner pair of yellow lines are fits to these calibrated uncertainties. The next sets of yellow lines show 2σ and 3σ uncertainties, scaled from the 1σ fit. Red squares highlight sources that are estimated to be ≥3σ variables. Bottom panel: same as the top panel, but now showing the variability significance metric, σ<sub>Δm, F140W</sub>, vs. m̄<sub>F140W</sub>. A linear scale is used for |σ<sub>Δm, F140W</sub>| ≤ 3 and a logarithmic scale elsewhere. The cross-matched known AGN in the HUDF (J. Lyu et al. 2022) are marked by open squares. Variable sources identified only by difference imaging are shown with open circles.</figcaption></figure>
<table><caption>**Table 2**
Number of Variable Sources Retrieved in Each Filter and the Combined Numbers</caption><thead><tr><th>Significance Level</th><th># F105W</th><th># F140W</th><th># F160W</th><th># Combined</th><th>f<sub>AGN</sub></th><th>F<sub>var</sub></th></tr></thead><tbody><tr><td>>2σ</td><td>222</td><td>191</td><td>185</td><td>521</td><td>15/31</td><td>2.07</td></tr><tr><td>>2.5σ</td><td>92</td><td>67</td><td>68</td><td>188</td><td>11/31</td><td>2.81</td></tr><tr><td>>3σ</td><td>43</td><td>34</td><td>39</td><td>99 + 10</td><td>10/31</td><td>3.10</td></tr></tbody></table>

Note. We also report the fraction of the 31 known AGN from the J. Lyu et al. (2022) sample, $f_{AGN}$, that are recovered. The inverse of this fraction is used as the variability completeness correction factor, $F_{var} = f_{AGN}^{-1}$.

JADES GTO program for spectroscopic redshifts using NIRSpec (A. J. Bunker et al. 2023), photometric redshifts from NIRCam + HST (M. J. Rieke et al. 2023), and the Ultraviolet Ultra Deep Field photometric catalog (M. Rafelski et al. 2015). Spectroscopic estimates, when available, are given preference over photometric ones. We note that the average photometric redshift uncertainty is of the order of ≲0.2, which justifies our choice of unit width bins when considering the overall redshift distribution (see next section).

### 3. Results

#### 3.1. Detected Variable Sources

Figure 1(a) shows the photometric variability ($\Delta m = m_1 - m_2$) of galactic nucleus sources in the HUDF in the F140W filter from 2012 (epoch 1) to 2023 (epoch 2), i.e., an 11 yr observer-frame time baseline, as a function of mean magnitude, $\bar{m}$, over the range from about 21 to 32 mag. A linear fit to the logarithm of the calibrated $1\sigma$ uncertainty estimates in 0.5 mag bins (yellow circles) is shown by the inner pair of yellow lines. Note that the yellow circles represent the mean photometric uncertainties as provided by SExtractor multiplied by an overall scale factor (0.84 for F140W) in order to match the inner 68th percentile of the photometric variation (black circles). This factor accounts for systematic uncertainties in the estimate of the photometric uncertainties, e.g., due to nonconstant point-spread function and correlated pixel noise. The fitted lines are then scaled by factors of two and three to show estimates of the $2\sigma$ and $3\sigma$ uncertainties. Sources that have values of $|\Delta m_{F140W}|$ that are greater than $3\sigma$ have been highlighted. Note that for the faintest sources with $m_{F140W} \gtrsim 30$ mag this corresponds to a variation of $\gtrsim 0.6$ mag.

Figure 1(b) shows the dimensionless variability metric $\sigma_{\Delta m, F140W}$ of each source versus mean magnitude. A linear scale is used for the inner range from -3 to +3 and then a logarithmic scaling at higher absolute values. The limits of $\sigma_{\Delta m, F140W} = 2, 2.5, 3$ are highlighted, showing the different levels of significance that we consider. From these data we identify 191, 67, and 39 sources that are >2$\sigma$, >2.5$\sigma$, and >3$\sigma$ variables, respectively, in the F140W observations (see Table 2).

In Figure 1(b) we also highlight the cross-match with the 31 previously known AGN from the sample of J. Lyu et al. (2022), with these sources having $\bar{m}_{F140W} \gtrsim 28$ mag. We see that four of these AGN have been detected as variable sources with $\sigma_{\Delta m} > 3$, still four with $\sigma_{\Delta m} > 2.5$, and five with $\sigma_{\Delta m} > 2$. In Section 3.3 we discuss our efficiency of recovering AGN via our global analysis of three measurements of variability via the F105W, F140W, and F160W observations.

The equivalent results of Figure 1 for the F105W and F160W observations are shown in the Appendix. From these data, the equivalent numbers of variable sources from the F105W (2009–2012, i.e., 3 yr time baseline) observation are

222, 92, and 43 sources at $2\sigma$, $2.5\sigma$, and $3\sigma$ variables, respectively. The F160W (2009–2012, i.e., 3 yr time baseline) observation yields 185, 68, and 39 sources at $2\sigma$, $2.5\sigma$, and $3\sigma$ variables, respectively. We see that the F105W, F140W, and F160W observations have similar efficiencies at detecting variable sources.

The difference imaging method identifies a further six variable sources, which were not detected as variable via the photometric method. In the following we consider these sources as part of the $\ge 3\sigma$ subsample. Additionally, three more sources previously identified at $2\sigma$ and one source identified at $2.5\sigma$ are promoted to $3\sigma$ for a total of 10 sources added at this significance level.

Combining the above results for all the filters, in the end we obtain a total sample of 521, 188, and 99 + 10 = 109 sources that show variability in at least one filter at the $2\sigma$, $2.5\sigma$, and $3\sigma$ levels, respectively. Of these, 482, 173, and 95 sources have been cross-matched with measured redshift catalogs, respectively, with the redshift distribution discussed in the next subsection. The photometric variability data for each source with a known redshift are presented in the electronic version of Table 2, with the displayed version here listing only the highest-redshift ($z > 6$) sources.

Finally, for completeness, we define the total variability significance of a source from the three observations using F105W, F140W, and F160W via

$$ \sigma_{\text{TOT}} = \sqrt{\sigma_{\text{F105W}}^2 + \sigma_{\text{F140W}}^2 + \sigma_{\text{F160W}}^2} \quad (2) $$

The values of $\sigma_{\text{TOT}}$ are also listed in Table 3. Note that, assuming independent Gaussian distributions for each filter measurement, the final $\sigma_{\text{TOT}}$ follows the statistics of a $\chi$ distribution. We checked for correlation between F105W and F160W variability, finding that they are largely uncorrelated. This may be caused by the visits in the 2009–2010 epochs being taken up to 1 yr apart (see Section 1). The values of $\sigma_{\text{TOT}}$ corresponding to significance levels of $2\sigma, 2.5\sigma$, and $3\sigma$ are 2.83, 3.30, and 3.76, respectively. Nevertheless, for definition of our samples of AGN candidates, we will use the metrics based on single filter variability significance, as described above.

#### 3.2. Redshift Distribution

Figure 2 presents a scatter plot of $\sigma_{\Delta m}$ versus $z$ for the F105W, F140W, and F160W data sets. We see that the redshift distribution of the sources in the sample extends out to $z \sim 9.5$, with the highest-redshift (at $2.5\sigma$) variability detections being at $z \sim 8$ (three sources). Among them, the highest-redshift variable source in our sample results at $z = 7.53^{+0.46}_{-0.29}$. Table 3 lists all 13 of the variable sources that we identify with $z > 6$ (note that the online electronic table contains the full sample at all redshifts).
<table><caption>**Table 3**

General properties of the $z > 6$ variable sources at $>2\sigma$ significance combining the three different filters. Values in bold text indicate whether the source shows a $\sigma$ threshold larger than 2 in a given filter.</caption><thead><tr><th>ID</th><th>R.A.<br>(deg)</th><th>Decl.<br>(deg)</th><th>m<sub>1,F105W</sub><br>(mag)</th><th>m<sub>2,F105W</sub><br>(mag)</th><th>m<sub>1,F160W</sub><br>(mag)</th><th>m<sub>2,F160W</sub><br>(mag)</th><th>m<sub>1,F140W</sub><br>(mag)</th><th>m<sub>2,F140W</sub><br>(mag)</th><th>σ<sub>F105W</sub><sup>c</sup></th><th>σ<sub>F160W</sub><sup>c</sup></th><th>σ<sub>F140W</sub><sup>d</sup></th><th>σ<sub>TOT</sub></th><th>Redshift</th></tr></thead><tbody><tr><td>2165</td><td>53.16823</td><td>-27.78223</td><td>30.00 ± 0.19</td><td>30.44 ± 0.17</td><td>30.53 ± 0.24</td><td>30.51 ± 0.35</td><td>30.74 ± 0.28</td><td>30.05 ± 0.16</td><td>-1.58</td><td>-0.02</td><td><strong>2.03</strong></td><td>2.57</td><td>7.53<sup>+0.46</sup><sub>-0.29</sub></td></tr><tr><td>3326</td><td>53.16744</td><td>-27.77195</td><td>30.08 ± 0.21</td><td>30.81 ± 0.24</td><td>30.86 ± 0.32</td><td>31.18 ± 0.64</td><td>30.88 ± 0.32</td><td>30.16 ± 0.17</td><td><strong>-2.31</strong></td><td>-0.49</td><td>1.85</td><td>3.0</td><td>7.27<sup>+0.11</sup><sub>-0.75</sub></td></tr><tr><td>2567</td><td>53.18626</td><td>-27.77897</td><td>29.31 ± 0.13</td><td>29.26 ± 0.06</td><td>28.79 ± 0.05</td><td>28.81 ± 0.07</td><td>28.70 ± 0.05</td><td>28.85 ± 0.06</td><td>0.49</td><td>-0.13</td><td><strong>-2.20</strong></td><td>2.26</td><td>7.27<sup>+0.11</sup><sub>-0.75</sub></td></tr><tr><td>1297</td><td>53.16481</td><td>-27.7882</td><td>29.18 ± 0.09</td><td>29.19 ± 0.05</td><td>28.71 ± 0.05</td><td>28.52 ± 0.06</td><td>28.70 ± 0.04</td><td>28.60 ± 0.04</td><td>0.16</td><td><strong>2.76</strong></td><td>1.52</td><td>3.16</td><td>7.227<sup>c</sup></td></tr><tr><td>3270</td><td>53.16553</td><td>-27.77259</td><td>29.57 ± 0.13</td><td>30.02 ± 0.12</td><td>29.60 ± 0.10</td><td>30.14 ± 0.24</td><td>29.67 ± 0.11</td><td>29.76 ± 0.12</td><td><strong>-2.34</strong></td><td><strong>-2.06</strong></td><td>-0.6</td><td>3.18</td><td>7.17<sup>+0.09</sup><sub>-0.07</sub></td></tr><tr><td>1511<sup>a</sup></td><td>53.16192</td><td>-27.78699</td><td>28.48 ± 0.06</td><td>28.22 ± 0.02</td><td>28.05 ± 0.02</td><td>27.77 ± 0.03</td><td>27.83 ± 0.02</td><td>27.99 ± 0.02</td><td><strong>4.58</strong></td><td><strong>7.41</strong></td><td><strong>-4.8</strong></td><td>9.95</td><td>6.74<sup>+0.04</sup><sub>-0.07</sub></td></tr><tr><td>2080</td><td>53.17732</td><td>-27.78238</td><td>28.21 ± 0.04</td><td>28.14 ± 0.02</td><td>28.04 ± 0.03</td><td>28.07 ± 0.04</td><td>28.06 ± 0.03</td><td>28.01 ± 0.03</td><td><strong>2.15</strong></td><td>-0.51</td><td>1.36</td><td>2.59</td><td>6.73<sup>+0.05</sup><sub>-0.03</sub></td></tr><tr><td>1058</td><td>53.15794</td><td>-27.79093</td><td>30.03 ± 0.19</td><td>30.73 ± 0.22</td><td>29.96 ± 0.14</td><td>30.28 ± 0.27</td><td>30.6 ± 0.15</td><td>30.14 ± 0.18</td><td><strong>-2.40</strong></td><td>-1.07</td><td>-0.36</td><td>2.65</td><td>6.56<sup>+0.001</sup><sub>-6.37</sub></td></tr><tr><td>101159<sup>b</sup></td><td>53.16052</td><td>-27.78593</td><td>28.96 ± 0.37</td><td>≥30.1</td><td>28.96 ± 0.37</td><td>≥30.1</td><td>≥29.9</td><td>≥29.9</td><td>...</td><td>...</td><td>...</td><td>...</td><td>6.54<sup>+2.45</sup><sub>-2.59</sub></td></tr><tr><td>305</td><td>53.1564</td><td>-27.80044</td><td>30.02 ± 0.19</td><td>29.96 ± 0.11</td><td>30.51 ± 0.24</td><td>29.87 ± 0.19</td><td>30.23 ± 0.18</td><td>30.17 ± 0.18</td><td>0.48</td><td><strong>2.08</strong></td><td>0.27</td><td>2.15</td><td>6.26<sup>+0.19</sup><sub>-0.02</sub></td></tr><tr><td>3391</td><td>53.16101</td><td>-27.77124</td><td>29.92 ± 0.18</td><td>29.76 ± 0.09</td><td>29.78 ± 0.12</td><td>30.50 ± 0.33</td><td>29.89 ± 0.13</td><td>30.00 ± 0.15</td><td>0.95</td><td><strong>-2.01</strong></td><td>-0.52</td><td>2.29</td><td>6.21<sup>+0.31</sup><sub>-0.05</sub></td></tr><tr><td>2736</td><td>53.18742</td><td>-27.77797</td><td>29.44 ± 0.15</td><td>30.53 ± 0.25</td><td>31.04 ± 0.50</td><td>31.72 ± 1.16</td><td>30.93 ± 0.46</td><td>30.37 ± 0.29</td><td><strong>-3.67</strong></td><td>0.0</td><td>0.97</td><td>3.79</td><td>6.08<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>2680</td><td>53.15549</td><td>-27.77824</td><td>29.97 ± 0.18</td><td>30.33 ± 0.15</td><td>31.43 ± 0.53</td><td>30.33 ± 0.28</td><td>30.09 ± 0.15</td><td>30.99 ± 0.37</td><td>-1.36</td><td>1.78</td><td><strong>-2.31</strong></td><td>3.22</td><td>6.01<sup>+0.05</sup><sub>-0.06</sub></td></tr></tbody></table>

Notes. <sup>a</sup> This source corresponds to source 1052123 in Paper I. <sup>b</sup> This source is taken from M. J. Hayes et al. (2024) and was only detected via the image subtraction method as described in the text. <sup>c</sup> Comparing epochs taken in 2008–9 vs. 2012. <sup>d</sup> Comparing epochs taken in 2012 vs. 2023.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Variability significance metric, $\sigma_{\Delta m}$, as a function of redshift for the three filters F105W, F140W, and F160W. The cross-matched known AGN in the HUDF (J. Lyu et al. 2022) are marked by open squares. Sources identified by difference imaging are shown with open circles.</figcaption></figure>

Figure 3(a) shows the redshift histogram of our sample of detected variable sources. Different colors depict the different samples according to the significance level (>2$\sigma_m$, 2.5$\sigma_m$, and 3$\sigma_m$) as detected in at least one filter. We see that most candidates reside at $z \lesssim 4$ for all cases of significance threshold. Figure 3(b) shows the fraction of galaxies that are detected to have variable nuclei. This fraction appears relatively constant with redshift and is seen to be elevated compared to the false-positive levels, with this enhancement relatively greater at higher values of $\sigma$. However, the absolute number of variable sources is often dominated by those found at 2$\sigma$, even after subtracting off the number of expected false positives (see the next subsection).

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top panel: redshift distribution of the variable candidate sources. Different colors refer to the specific significance thresholds for variability detection in units of $\sigma$. The histogram includes sources selected both via variability and via difference imaging in the three filters. Bottom panel: fraction of galaxies detected as variable sources as a function of the redshift at different $\sigma$ levels (horizontal bars). The width of the bars indicates the adopted binning in redshift. Vertical error bars depict Poisson uncertainties, while dotted lines show the expected number of false positives given assuming a combined Gaussian statistic in the three filters (see text).</figcaption></figure>

For a visual representation of the high-redshift variable sources in the field, Figures 4 and 5 show their positions superimposed on the HUDF image and cutout frames in different bands used for the visual inspection, respectively. Different symbols indicate the different ranges of significance level in units of $\sigma$ as detailed in the legend, color-coded according to the estimated redshift. Zoom-in images of the highest-redshift ($z > 7$) variable candidates are presented in small inset panels in Figure 4.

### 3.3. Comoving Number Density of SMBHs

Next, assuming that the detected variable sources are AGN powered by accreting SMBHs, we use the numbers of variable sources to estimate the comoving number density of SMBHs, $n_{SMBH}(z)$. As discussed in Paper I, we expect variable sources detected at high redshifts to be more likely to be AGN rather than SNe or other stellar transients. In addition, our focus on galactic nuclei will also tend to reduce SN interlopers in our AGN sample. Very high redshift sources, i.e., $z \gtrsim 6$, where the lower $1\sigma$ uncertainty in redshift is compatible with a low-redshift interloper scenario, i.e., $z_{phot} - z_{low} \gtrsim 2$, are discarded from the number density estimate.

Figure 6 and Table 4 present the information related to this analysis, starting with the adopted redshift intervals and the
<figure><img src="image_6.png" /><figcaption>Figure 4. Angular distribution of variable candidates in the HUDF at $z \ge 6$. The estimated redshift is shown by the color of the symbol with reference to the color bar. Different symbols refer to different significance levels (see legend). We also highlight sources at $z > 7$ with zoom-in thumbnails.</figcaption></figure>

raw number of SMBHs in each interval, $N_{SMBH,raw}$ (for each of the $2\sigma$, $2.5\sigma$, and $3\sigma$ significance levels). Next, we list the total number of galaxies analyzed, $N_{gal}$, which is used to estimate the number of false positives, $N_{fp}$, that would result at each significance level. Note that there are three chances for selection, so that the fraction of false positives $f_{fp}$, assuming uncorrelated Gaussian distributions in the three filters, is given by

$$f_{fp} = 1 - (1 - f_{fp}, \sigma)^3, \qquad (3)$$

where $f_{fp,\sigma}$ is the fraction of false-positive detections expected according to a Gaussian distribution at a given $\sigma$ level. Evaluating this equation at the $2\sigma$, $2.5\sigma$, and $3\sigma$ levels, we find $f_{fp}$ values of 0.13, 0.037, and 0.008, respectively. The final estimate of the number of SMBHs in each interval is then $N_{SMBH,raw-fp} = N_{SMBH,raw} - N_{fp}$. Then, using the comoving volume of the HUDF footprint projected across each redshift interval, $V$, we evaluate the comoving number densities $n_{SMBH,raw}$ and $n_{SMBH,raw-fp}$.

We consider $n_{SMBH,raw-fp}$ to be the most reliable, direct estimate of a lower limit to the comoving number density of SMBHs. These values are plotted as open squares in Figure 6. We see that in the lowest-redshift bin from $z = 0$ to $0.5 n_{SMBH,raw-fp} \simeq 6 \times 10^{-3} \text{cMpc}^{-3}$ regardless of the $\sigma$ threshold. We note that these values are a few times higher than the previous $z \sim 0$ estimate of $n_{SMBH}$ estimated by M. Vika et al. (2009) and N. Banik et al. (2019). At the highest redshifts, i.e., $z \sim 7-8$, we find $n_{SMBH,raw-fp} \sim 10^{-4} \text{cMpc}^{-3}$. In between these redshift extremes we derive smoothly declining intermediate values of $n_{SMBH,raw-fp}$ as redshift increases.

However, there are at least two factors causing the estimate of $n_{SMBH,raw-fp}$ to be an incomplete census of SMBHs. First, we do not expect all AGN to vary during a given time interval, especially if that interval becomes relatively short in the AGN rest frame. To estimate a variability incompleteness factor, $F_{var}$, we compare our sources with the sample of 31 known AGN in the GOODS-S field, as compiled by J. Lyu et al. (2022). These AGN were identified through various methods, including mid-IR colors (four AGN), X-ray luminosities (seven AGN), radio-loudness (seven AGN), optical spectroscopy (one
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Cutout images of the sources listed in Table 3. In particular for each source, from left to right, we have an F105W difference image, F140W difference image, and median ACS image (F606W + F775W + F814W for $z > 6.9$ or F435W + F606W + F775W for $z < 6.9$) for identifying dropouts/ interlopers; the median F105W image (median stack of both epochs); the median F140W image; and the JADES F277W image (deepest JADES filter image).</figcaption></figure>

AGN), and, in some cases, variability. Note that several AGN were identified by more than one diagnostic technique. Table 2 reports the recovery fraction, $f_{AGN}$, of these 31 known AGN in our 2$\sigma$, 2.5$\sigma$, and 3$\sigma$ significance samples, i.e., 15/31, 11/31, and 10/31, respectively (note that the 3$\sigma$ result includes the difference-image-selected sources). As a simple empirical method, we adopt the inverse of these recovery fractions as our variability incompleteness factors, $F_{var}$, i.e., taking values of

2.07, 2.81, and 3.10. We note that these factors have significant Poisson uncertainties, which we include in the analysis. We also note that the J. Lyu et al. (2022) sample is biased to lower redshifts, i.e., $z \lesssim 3$. We expect that variability incompleteness would tend to become larger at higher redshifts owing to shorter rest-frame time baselines. On the other hand, the variability amplitude is expected to grow progressively at bluer wavelengths, i.e. those that are observed at higher redshifts.

Our approach does not account for sources of incompleteness arising from the diverse AGN selection methods employed to compile the J. Lyu et al. (2022) sample. For instance, dust-obscured AGN might remain undetected in X-ray surveys, or inadequate spectral resolution could hinder precise AGN diagnostics. Additionally, host galaxy contamination, which varies with redshift and host galaxy luminosity, could influence the number of variable sources recovered via variability studies, as well as the observed inverse relationship between variability significance and AGN luminosity (e.g., K. L. Suberlak et al. 2021). Nevertheless, we argue that incorporating these complexities would statistically increase the expected number of recovered AGN in the field, of which only a fraction would be varying in the optical, consequently leading to lower recovered fractions (i.e., larger incompleteness corrections). In this scenario, our $F_{var}$ would thus represent, despite considerable uncertainties, a conservative estimate. More extensive, longer-term monitoring of AGN variability across a range of luminosities, wavelengths, and redshifts is needed for improved estimates of $F_{var}$.

After the above correction for variability incompleteness, in the low-$z$ interval we obtain comoving number densities of SMBHs of $n_{SMBH, var} \simeq 2 \times 10^{-2}$ cMpc$^{-3}$ (open diamonds in Figure 6), with little sensitivity to the choice of significance level. This may indicate that our estimate of the variability incompleteness factor as a function of significance level is reasonable. The smallest uncertainties result from the 2.5$\sigma$-selected sample, formally yielding $n_{SMBH, var} = (1.7 \pm 1.4) \times 10^{-2}$ cMpc$^{-3}$. We regard this estimate as our best measure of the true local number density of SMBHs. However, it is likely to still be a lower limit if there are very faint AGN whose flux and flux variations would not be detected by the HST observations. At the highest redshifts, $z = 7-8$, we find $n_{SMBH, var} \simeq 2 \times 10^{-4}$ cMpc$^{-3}$. Given that we have adopted a redshift-independent variability correction factor, the redshift dependence of $n_{SMBH, var}$ is the same as that of $n_{SMBH, raw-fp}$.

Finally, we correct for “luminosity incompleteness” via a luminosity correction factor, $F_{lum}$, i.e., $n_{SMBH, lum} = F_{lum}n_{SMBH, var}$. In particular, at high redshift we are only able to detect relatively bright AGN (e.g., at $z \sim 6-7$ this corresponds to $M_{UV} \simeq -18.6$), so this correction factor will boost the inferred number density of SMBHs. Following the method of Paper I and Y. Harikane et al. (2023a), who detected AGN via broad emission lines, we extrapolate from our detected AGN down an assumed luminosity function to a level of $M_{UV} = -17$ mag. We note that this is an arbitrary level, but by adopting this choice we are able to make a fair comparison to these previous studies. The value of the luminosity correction factor, $F_{lum}$, generally becomes larger with redshift, but it also depends on the form of the UV luminosity function (UVLF) we assume. Specifically, we integrate the double power-law fits detailed in S. L. Finkelstein & M. B. Bagley (2022) from redshift $z \sim 3$ to 9. At lower redshift, we integrate Schechter function fits to the UVLF from UVCANDELS and GALEX survey programs (S. Arnouts et al. 2005; L. Sun et al. 2024). At
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Evolution of the comoving number density of SMBHs, $n_{SMBH}$, vs. redshift, $z$, in intervals marked by the vertical dotted lines. We report constraints based on candidate AGN detected as $2\sigma$, $2.5\sigma$, and $3\sigma$ variability significance thresholds with turquoise, green, and magenta symbols, respectively. In each bin, we report number densities based on raw counts (crosses), raw counts corrected for expected false positives (open squares), variability-incompleteness-corrected counts (open diamonds), and luminosity-incompleteness-corrected counts (filled diamonds; see text). Population III.1 SMBH seeding models (N. Banik et al. 2019; J. Singh et al. 2023) with isolation distance parameters, $d_{iso}$, from 50 to 200 kpc (proper distance) are shown by the colored solid lines. At low redshifts these decrease compared to the maximum value attained (dotted lines) owing to mergers. The green dashed line shows the SMBH seeding assumed in M. Vogelsberger et al. (2014) based on an HMT above $7.1 \times 10^{10} M_{\odot}$ (shaded region shows a factor of two variation in this mass scale). The SMBH abundance achieved in a simulation of SMBH seeding via DC (S. Chon et al. 2016) is shown by the pink dashed line. The dashed light-blue line shows the black hole number density as obtained in the Romulus25 simulation based on the physical condition of baryon particles (M. Tremmel et al. 2017).</figcaption></figure>

$z \lesssim 4$ we see that $F_{\text{lum}}$ becomes smaller than unity as our detection limit becomes fainter than the UVLF lower bound of $M_{\text{UV}} = -17$. Note that in the low-$z$ interval, i.e., at $z = 0-0.5$, our sensitivity corresponds to an absolute magnitude of about $M_{\text{UV}} = -8.5$ assuming an average redshift of $z \sim 0.25$. This explains why we consider $n_{\text{SMBH, var}}$ as our best estimate at low redshifts. In Table 4 we report the values of $F_{\text{lum}}$ used in each redshift interval.

After the above correction factor for luminosity incompleteness, in the high-redshift regime ($z = 7-8$) we find $n_{\text{SMBH, lum}} \sim (3 \pm 6) \times 10^{-3} \text{cMpc}^{-3}$. We caution that this result is consistent with a null density and involves a relatively large luminosity function correction factor of $F_{\text{lum}} \simeq 13$ and that the systematic uncertainty associated with this correction factor is not included in the above estimate. The redshift 6–7 interval has similarly high values of $n_{\text{SMBH, lum}}$ (although the sample based on $2\sigma$ significance does not yield a number of sources greater than those expected from false positives). At intermediate redshifts, $z = 3-6$, we notice that $n_{\text{SMBH, lum}}$ takes smaller values $\sim 10^{-3}$ to $10^{-4} \text{cMpc}^{-3}$. We do not expect such a rapid decline in the true number density of SMBHs, so this result may indicate either that a larger fraction of SMBHs were brighter than $M_{\text{UV}} = -17$ mag at $z > 6$ or that the values at high redshift are somewhat overestimated via their luminosity function correction factor. At low redshifts, we see that $n_{\text{SMBH, lum}}$ rises again to values of

$\lesssim 10^{-3} \text{cMpc}^{-3}$. However, at $z = 0-0.5$ this is a factor of about 20 smaller than the number we infer after variability completeness correction. This is simply a reflection of the fact that the observations are able to detect AGN that are much fainter than the arbitrary limit of $M_{\text{UV}} = -17$ mag.

In Figure 7 we show our $3\sigma$ results for $n_{\text{SMBH, var}}$ and $n_{\text{SMBH, lum}}$, along with several other estimates of $n_{\text{SMBH}}$. As mentioned, the previous estimates of $n_{\text{SMBH}}(z=0)$ of $\sim 5 \times 10^{-3} \text{cMpc}^{-3}$ (N. Banik et al. 2019) and $\sim 9 \times 10^{-3} \text{cMpc}^{-3}$ (M. Vika et al. 2009) now appear to be superseded by about a factor of a few with our HUDF variability study estimate of $n_{\text{SMBH, var}}$. Similarly, we see that our values of $n_{\text{SMBH, var}}$ out to $z = 3$ are about 2 times higher than the values of $n_{\text{SMBH}}$ implied by the 31 previously known AGN in the HUDF (J. Lyu et al. 2022).

At high redshifts we find consistency with our result from Paper I based on three detected AGN candidates between $z = 6$ and 7. Our results for $n_{\text{SMBH, lum}}$ at $z > 6$ exceed those inferred from the JWST-selected sample of Y. Harikane et al. (2023a) by about a factor of 5. From $z = 4$ to 6 our estimates of $n_{\text{SMBH, lum}}$ are instead about 5 times lower than those of Y. Harikane et al. (2023a). For all these high-$z$ results, there are large uncertainties due to the small number of directly detected sources. Indeed, our $2\sigma$-selected sample at $z > 7$ has the smallest Poisson errors.
<table><caption>**Table 4**

Census of SMBHs across the Universe in Various Redshift Intervals Out to z = 9</caption><thead><tr><th>Redshift</th><th>σ</th><th>N<sub>SMBH,raw</sub></th><th>N<sub>gal</sub></th><th>N<sub>fp</sub></th><th>N<sub>SMBH,raw-fp</sub><br>(10<sup>3</sup> cMpc<sup>-3</sup>)</th><th>V<br>(10<sup>-3</sup> cMpc<sup>-3</sup>)</th><th>n<sub>SMBH,raw</sub><br>(10<sup>-3</sup> cMpc<sup>-3</sup>)</th><th>n<sub>SMBH,raw-fp</sub><br>(10<sup>-3</sup> cMpc<sup>-3</sup>)</th><th>δn<sub>SMBH,raw-fp</sub><br>(10<sup>-3</sup> cMpc<sup>-3</sup>)</th><th>n<sub>SMBH,var</sub><br>(10<sup>-3</sup> cMpc<sup>-3</sup>)</th><th>δn<sub>SMBH,var</sub></th><th>F<sub>lum</sub><br>(10<sup>-3</sup> cMpc<sup>-3</sup>)</th><th>n<sub>SMBH,lum</sub><br>(10<sup>-3</sup> cMpc<sup>-3</sup>)</th><th>δn<sub>SMBH,lum</sub></th></tr></thead><tbody><tr><td>0-0.5</td><td>2</td><td>55</td><td>...</td><td>52.1</td><td>27.3</td><td>...</td><td>58.5</td><td>6.4</td><td>7.7</td><td>13.2</td><td>16.2</td><td>...</td><td>0.6</td><td>0.7</td></tr><tr><td>...</td><td>2.5</td><td>21</td><td>400</td><td>14.7</td><td>16.1</td><td>1.0</td><td>20.8</td><td>6.1</td><td>4.8</td><td>17.2</td><td>14.5</td><td>0.046</td><td>0.8</td><td>0.7</td></tr><tr><td>...</td><td>3</td><td>10</td><td>...</td><td>3.2</td><td>10.7</td><td>...</td><td>9.9</td><td>6.7</td><td>4.9</td><td>20.7</td><td>16.4</td><td>...</td><td>1.0</td><td>0.8</td></tr><tr><td>0.5-1</td><td>2</td><td>86</td><td>...</td><td>66.5</td><td>29.9</td><td>...</td><td>19.9</td><td>4.4</td><td>2.2</td><td>9.0</td><td>5.1</td><td>...</td><td>0.7</td><td>0.4</td></tr><tr><td>...</td><td>2.5</td><td>32</td><td>510</td><td>18.7</td><td>24.0</td><td>4.3</td><td>7.4</td><td>3.0</td><td>1.5</td><td>8.5</td><td>4.9</td><td>0.074</td><td>0.6</td><td>0.4</td></tr><tr><td>...</td><td>3</td><td>23</td><td>...</td><td>4.1</td><td>14.8</td><td>...</td><td>5.3</td><td>4.4</td><td>2.4</td><td>13.5</td><td>8.6</td><td>...</td><td>1.0</td><td>0.6</td></tr><tr><td>1 2</td><td>2</td><td>162</td><td>...</td><td>118.2</td><td>50.5</td><td>...</td><td>11.0</td><td>3.0</td><td>0.9</td><td>6.1</td><td>2.4</td><td>...</td><td>1.5</td><td>0.4</td></tr><tr><td>...</td><td>2.5</td><td>70</td><td>907</td><td>33.4</td><td>32.5</td><td>14.9</td><td>4.7</td><td>2.5</td><td>0.7</td><td>6.9</td><td>2.9</td><td>0.148</td><td>1.7</td><td>0.4</td></tr><tr><td>...</td><td>3</td><td>38</td><td>...</td><td>7.3</td><td>18.7</td><td>...</td><td>2.4</td><td>1.9</td><td>0.8</td><td>6.0</td><td>3.2</td><td>...</td><td>1.5</td><td>0.5</td></tr><tr><td>2-3</td><td>2</td><td>92</td><td>...</td><td>82.4</td><td>-0.5</td><td>...</td><td>5.5</td><td>0.6</td><td>1.1</td><td>...</td><td>1.2</td><td>...</td><td>0.6</td><td>0.5</td></tr><tr><td>...</td><td>2.5</td><td>26</td><td>632</td><td>23.3</td><td>7.7</td><td>16.8</td><td>1.5</td><td>0.2</td><td>0.3</td><td>0.4</td><td>0.9</td><td>0.405</td><td>0.3</td><td>0.4</td></tr><tr><td>...</td><td>3</td><td>13</td><td>...</td><td>5.2</td><td>11.9</td><td>...</td><td>0.8</td><td>0.5</td><td>0.3</td><td>1.5</td><td>1.0</td><td>...</td><td>0.8</td><td>0.4</td></tr><tr><td>3-4</td><td>2</td><td>37</td><td>...</td><td>40.5</td><td>7.5</td><td>...</td><td>2.3</td><td>-0.2</td><td>0.4</td><td>-0.3</td><td>0.8</td><td>...</td><td>-0.2</td><td>0.6</td></tr><tr><td>...</td><td>2.5</td><td>12</td><td>303</td><td>11.2</td><td>4.8</td><td>15.9</td><td>0.8</td><td>0.1</td><td>0.2</td><td>0.2</td><td>0.6</td><td>0.72</td><td>0.1</td><td>0.4</td></tr><tr><td>...</td><td>3</td><td>2</td><td>...</td><td>2.4</td><td>1.6</td><td>...</td><td>0.1</td><td>0.0</td><td>0.1</td><td>-0.1</td><td>0.3</td><td>...</td><td>-0.1</td><td>0.2</td></tr><tr><td>4-5</td><td>2</td><td>21</td><td>...</td><td>22.8</td><td>1.2</td><td>...</td><td>1.5</td><td>-0.1</td><td>0.3</td><td>-0.3</td><td>0.7</td><td>...</td><td>-0.4</td><td>0.9</td></tr><tr><td>...</td><td>2.5</td><td>5</td><td>175</td><td>6.4</td><td>1.6</td><td>14.4</td><td>0.3</td><td>-0.1</td><td>0.2</td><td>-0.3</td><td>0.5</td><td>1.37</td><td>-0.4</td><td>0.6</td></tr><tr><td>...</td><td>3</td><td>4</td><td>...</td><td>1.4</td><td>2.6</td><td>...</td><td>0.3</td><td>0.2</td><td>0.2</td><td>0.6</td><td>0.7</td><td>...</td><td>0.8</td><td>0.9</td></tr><tr><td>5-6</td><td>2</td><td>16</td><td>...</td><td>15.6</td><td>3.4</td><td>...</td><td>1.2</td><td>0.0</td><td>0.3</td><td>0.1</td><td>0.6</td><td>...</td><td>0.2</td><td>1.8</td></tr><tr><td>...</td><td>2.5</td><td>2</td><td>120</td><td>4.4</td><td>2.6</td><td>12.9</td><td>0.2</td><td>-0.2</td><td>0.1</td><td>-0.5</td><td>0.4</td><td>2.8</td><td>-1.5</td><td>1.2</td></tr><tr><td>...</td><td>3</td><td>1</td><td>...</td><td>1.0</td><td>1.0</td><td>...</td><td>0.1</td><td>0.0</td><td>0.1</td><td>0.0</td><td>0.2</td><td>...</td><td>0.0</td><td>0.7</td></tr><tr><td>6-7</td><td>2</td><td>8</td><td>...</td><td>13.9</td><td>-1.8</td><td>...</td><td>0.7</td><td>-0.5</td><td>0.3</td><td>-1.1</td><td>0.6</td><td>...</td><td>-6.2</td><td>3.7</td></tr><tr><td>...</td><td>2.5</td><td>4</td><td>107</td><td>4.0</td><td>1.9</td><td>11.5</td><td>0.3</td><td>0.0</td><td>0.2</td><td>0.0</td><td>0.5</td><td>5.8</td><td>0.1</td><td>2.8</td></tr><tr><td>...</td><td>3</td><td>4</td><td>...</td><td>0.9</td><td>4.1</td><td>...</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.8</td><td>1.1</td><td>...</td><td>4.9</td><td>6.3</td></tr><tr><td>7-8</td><td>2</td><td>5</td><td>...</td><td>3.9</td><td>6.7</td><td>...</td><td>0.5</td><td>0.1</td><td>0.2</td><td>0.2</td><td>0.5</td><td>...</td><td>2.8</td><td>6.0</td></tr><tr><td>...</td><td>2.5</td><td>1</td><td>30</td><td>1.1</td><td>2.2</td><td>10.9</td><td>0.1</td><td>0.0</td><td>0.1</td><td>0.0</td><td>0.3</td><td>12.9</td><td>-0.4</td><td>3.5</td></tr><tr><td>...</td><td>3</td><td>0</td><td>...</td><td>0.3</td><td>1.6</td><td>...</td><td>0.0</td><td>0.0</td><td>0.1</td><td>-0.1</td><td>0.2</td><td>...</td><td>-1.2</td><td>2.2</td></tr></tbody></table>

**Note.** In each interval we report the raw number of variable sources, $N_{SMBH,raw}$, at $2\sigma$, $2.5\sigma$, and $3\sigma$ significance (see text), based on analysis of a total number of galaxies, $N_{gal}$, in that interval. Next, we list the number of expected false positives, $N_{fp}$. The excess number of variables above this false-positive level is listed as $N_{SMBH,raw-fp}$. The comoving volume for each interval is listed next, which is then used to evaluate comoving number densities based on raw ($n_{SMBH,raw}$) and false-positive-corrected ($n_{SMBH,raw-fp}$) counts. For the latter we also list its uncertainty due to Poisson counting statistics, $\delta n_{SMBH,raw-fp}$. Next, we list the variability-incompleteness-corrected estimate $n_{SMBH,var}$ and its uncertainty (see text). Next, we list the luminosity incompleteness correction factor, $F_{lum}$ (see text), which is then used to estimate $n_{SMBH,lum}$ and its uncertainty.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Same as Figure 6, but now only showing $n_{\text{SMBH}}$ estimated from $3\sigma$ variables. Various other observational constraints are displayed. The red diamond shows the observational constraint derived in Paper I (M. J. Hayes et al. 2024) based on three detected AGN between $z=6$ and 7. Broad emission line sources at $z=4-7$ (Y. Harikane et al. 2023a) are shown with the black diamonds, which use similar luminosity correction factors to our study. Brown stars indicate the number densities from the sample of 31 previously known AGN in the HUDF (J. Lyu et al. 2022). The black square shows a previous estimate of the local ($z=0$) comoving number density of SMBHs (N. Banik et al. 2019).</figcaption></figure>

### 3.4. Implications for SMBH Seeding Models

Our estimates of $n_{\text{SMBH}}$ shown in Figures 6 and 7 place strong constraints on SMBH seeding models. At high redshifts, $z>6$, we see that our estimates of $n_{\text{SMBH,lum}} \gtrsim 10^{-3} \text{cMpc}^{-2}$ are about a factor of 10 higher than the fiducial HMT model used in the Illustris family of simulations (M. Vogelsberger et al. 2014), in which dark matter halos with masses $>7 \times 10^{10} \text{M}_\odot$ are seeded with SMBHs. The range of the green band shown in Figures 6 and 7 corresponds to varying this threshold mass by factors of two. Thus, an HMT model with a threshold mass of $\sim 3 \times 10^{10} \text{M}_\odot$, which would be consistent with the $z=0$ constraints, would be closer to the high-$z$ data, but still low by about a factor of few.

Figures 6 and 7 also compare our constraints with the Romulus25 simulation (M. Tremmel et al. 2017). We see that this simulation, which seeds SMBHs using simple threshold conditions of gas properties (see Section 1), produces a high number density of SMBHs at early times, i.e., $n_{\text{SMBH}} \sim 5 \times 10^{-2} \text{cMpc}^{-3}$ at $z \sim 12$. This rises by almost an order of magnitude by $z \simeq 4$, before then undergoing a modest decline by $z=0$. Given its high number densities of SMBHs, this model remains consistent with our observational constraints on $n_{\text{SMBH}}$.

As discussed in Section 1, DC seeding models struggle to reproduce the overall number densities of SMBHs. The most

optimistic of the DC models discussed in Section 1 is that of S. Chon et al. (2016) with $n_{\text{SMBH}} \sim 10^{-4} \text{Mpc}^{-3}$. As shown in Figures 6 and 7, this is about a factor of 50 below $n_{\text{SMBH}}$ inferred from our HUDF variability study. Furthermore, the later simulations of DC by J. H. Wise et al. (2019) and M. A. Latif et al. (2022) find values of $n_{\text{SMBH}} \lesssim 10^{-6} \text{Mpc}^{-3}$, i.e., a factor of at least $10^3$ below our high-$z$ estimate. The constraints are even more severe when considering our low-$z$ estimates of $n_{\text{SMBH}}$, which have smaller uncertainties and values of $\sim 2 \times 10^{-2} \text{Mpc}^{-3}$.

A striking feature of our low-$z$ and high-$z$ estimates of $n_{\text{SMBH}}$ is the relative similarity of $n_{\text{SMBH,lum}}$ at high $z$ with $n_{\text{SMBH, var}}$ at low $z$. There are a number of reasons to expect that $n_{\text{SMBH,lum}}(z>6)$ is a lower limit, i.e., since it only considers AGN down to $M_{\text{UV}} = -17$ mag, and that variability incompleteness factors may be larger given the shorter rest-frame time intervals that are probed by the HUDF observations. This would then indicate a relatively constant evolution in $n_{\text{SMBH}}$ from $z \sim 8$ to $z \sim 0$. Such an evolution is a key feature and prediction of the Population III.1 SMBH seeding model of N. Banik et al. (2019), which has been further explored by J. Singh et al. (2023) and V. Cammelli et al. (2025). The main parameter of the Population III.1 model is the isolation distance, $d_{\text{iso}}$, needed for minihalos to be Population III.1 sources that seed SMBHs (as opposed to Population III.2 sources that do not). Expressed in proper
distances at the time of formation, our $n_{SMBH}$ results from the HUDF favor models with $d_{iso} \lesssim 50$ kpc. This would correspond to comoving separations of $\lesssim 1$ Mpc at typical formation redshifts of $z \sim 20$. We note that this result is consistent with that of V. Cammelli et al. (2025) of $d_{iso} \lesssim 75$ kpc based on the galaxy stellar mass function and the SMBH occupation fraction as a function of the stellar mass.

As discussed by J. C. Tan et al. (2024), more realistic physical models for $d_{iso}$ based on R-type H II region expansion around supermassive Population III.1 stars that are expected to be the direct progenitors of the SMBHs yield estimates of the comoving number density of SMBHs of

$$n_{SMBH} = \frac{3}{4\pi R_{R}^{3}} \rightarrow 0.18 t_{*,10}^{-1} S_{53}^{-1} \text{ cMpc}^{-3}, \quad (4)$$

where $t_{*,10} = t_{*}/10\text{Myr}$, $t_{*}$ is the lifetime of the supermassive star, $S_{53} \equiv S/10^{53} \text{s}^{-1}$, $S$ is the H-ionizing luminosity of the star, and $R_{R}$ is the radius of R-type H II region expansion achieved in the intergalactic medium. This estimate can be viewed as an upper limit since it has assumed maximal close packing of the sources and ignored other contributions to the ionizing background from Population III.2, Population II, and AGN sources. It is interesting that our estimate of $n_{SMBH}$ at low $z$ from the HUDF is within a factor of a few of this fiducial estimate.

Similar to the analysis carried out in Paper I, a valuable test for SMBH seeding schemes consists in measuring the actual distance between candidate AGN pairs. It is helpful to remind here that the Population III.1 model with $d_{iso} \simeq 100$ kpc in proper distance would correspond to a comoving separation of about 3 cMpc (if the relative motion by $z \sim 7$ is negligible). In addition, the HUDF footprint at this epoch is 5.5 cMpc on a side. In our sample, at $z > 6$ we are left with the 25 sources listed in Table 3, and three of them have spectroscopic redshifts available. Since we cannot assess the true distance with even one photometric redshift in a single pair (already a $\Delta_{z} = 0.01$ at $z \sim 7$ would be $\sim 2.5$ cMpc, comparable to $d_{iso}$), we have searched among the pairs having both spectroscopic redshift available. In this high-$z$ regime, we find none, and so no constraint can be used to rule out specific scenarios. Obtaining spectroscopic redshifts of all the high-$z$ SMBH candidates is needed to achieve better constraints.

### 3.5. Comparison with Other AGN Diagnostics

Despite recent advances, it remains immensely challenging to obtain an accurate census of AGN at $z \gtrsim 6$ at faint magnitudes. Hot dust tracers (e.g., D. Stern et al. 2005) are shifted to the far-IR, and X-ray and radio facilities are not sufficiently sensitive to identify individual faint AGN at these distances. High-ionization UV emission lines (e.g., He II $\lambda$1640, C IV $\lambda$1550, N V $\lambda$1240) are sometimes observed in luminous targets (e.g., N. Laporte et al. 2017; R. Mainali et al. 2017), but they can be inconclusive, as these lines also form in low-metallicity starbursts (D. A. Berg et al. 2019; A. Saxena et al. 2020; P. Senchyna et al. 2020). JWST-NIRSpec microshutter array (MSA) observations have delivered deep IR spectroscopy, and several studies (e.g., A. J. Cameron et al. 2023; K. Boyett et al. 2024; W. Hu et al. 2024) could identify AGN via line flux ratio classification (BPT-like diagrams; J. A. Baldwin et al. 1981), although the targeted approach of

MSA spectroscopy is subject to photometric preselection biases and tends to be focused more toward the brighter end of the luminosity function. In addition, NIRSpec surveys face challenges in achieving completeness at $M_{UV} \simeq -17$ to $-20$, where galaxies are $\sim 100$ times more abundant.

A number of studies (e.g., S. H. Cohen et al. 2006; E. Pouliasis et al. 2019; R. O'Brien et al. 2024) have shown the effectiveness of using photometric monitoring with HST to identify variable sources, finding large numbers of AGN at intermediate luminosities, that evade other selection techniques (X-ray, radio, IR colors). In addition, as demonstrated in this series of papers, we advocate that time variability is an efficient observable for AGN identification in the high-redshift Universe.

### 3.6. Supernova Contamination

Our observations are sensitive to the presence of stellar transients, especially SNe, which may be a source of contamination in our counting of AGN. In Paper I we identified three variable sources that are SNe, two of which are apparently hostless with indeterminable redshifts. Supernovae are expected to follow either the stellar mass distribution (thermonuclear) or star formation rate (SFR) distribution (core collapse). Statistically in a galaxy these are centrally concentrated, so our approach of focusing on deep-stack centroiding in galactic nuclei, while helping to avoid off-nuclear SNe, would still be subject to some level of stellar transient contaminants. We note that we do not consider tidal disruption events, i.e., stellar disruption by an SMBH, to be contaminants in our census.

As discussed in Paper I, our HST images are not deep enough to detect ordinary core-collapse SNe at $z \gtrsim 3$. We are potentially sensitive to thermonuclear SNe at these redshifts if they exist, but the delay times for the main Type Ia channels mean that they should be very rare at high $z$. Superluminous SNe are an interesting possibility and could be detected out to $z \sim 6$. However, their rates should follow the cosmic SFR density, and we expect very low numbers in the small volume probed by the HUDF. A follow-up paper in this series will focus on the full census of detected SN candidates.

## 4. Conclusions

Through photometric monitoring conducted in three epochs (2008/2009, 2012, and 2023) we have discovered many variable sources in the HUDF, which we consider to be likely AGN candidates tracing the presence of SMBHs. We advocate that variability searches are a highly effective and comprehensive tool for identifying AGN in deep imaging surveys. While only a fraction of AGN can be detected through their variability, the significant advantage of this approach lies in its ability to survey the entire field with no prior selection. Unlike other methods such as radio or X-ray diagnostics, which are limited by detection thresholds, variability searches ensure that any object captured in imaging can be tested for its activity.

Historically, high-redshift luminous quasars have been known for some time (X. Fan et al. 2006; D. J. Mortlock et al. 2011; E. Bañados et al. 2018), though they are typically much brighter—by four magnitudes or more—than the sources we detect. More recently, there have been reports of AGN at even higher redshifts with luminosities closer to $L^{*}$ (e.g., Y. Harikane et al. 2023a; R. L. Larson et al. 2023;
R. Maiolino et al. 2024). However, our survey, which was restricted to the small volume of the HUDF, has the ability to place stronger constraints than any previous work on the number densities of SMBHs at these redshifts. We have shown how the 13 variable sources we detect between $z = 6$ and 8 imply high comoving number densities of $n_{SMBH} \gtrsim 2 \times 10^{-3} \text{cMpc}^{-3}$. In addition, our results at low redshift make a new measurement of $n_{SMBH} \sim 2 \times 10^{-2} \text{cMpc}^{-3}$ in this regime that is about a factor of 5 higher than previous estimates.

Our measurements constrain SMBH seeding mechanisms and their implementation in cosmological simulations. For example, our estimate of $n_{SMBH}$ at $z = 7-8$ is about a factor of 10 higher than expected in the fiducial HMT model used in the Illustris family of simulations (M. Vogelsberger et al. 2014). Our estimates, including at $z \sim 0$, are $\sim 10^2-10^4$ times larger than produced in simulations of SMBH seeding via DC (S. Chon et al. 2016; J. H. Wise et al. 2019; M. A. Latif et al. 2022).

On the other hand, the relative constancy of $n_{SMBH}$ from high to low $z$ confirms a prediction of the Population III.1 SMBH seeding model of N. Banik et al. (2019) and constrains its main isolation distance parameter, $d_{iso}$, to be $\lesssim 75$ kpc. Such values are also consistent with recent estimates of SMBH occupation fractions as a function of the stellar mass (V. Cammelli et al. 2025). They are also consistent with fiducial Population III.1 model expectations of ionizing feedback for the isolation distance (J. C. Tan et al. 2024).

Future papers in this series will examine the properties of the host galaxies of our detected AGN candidates (A.R. Young et al. 2025) and carry out a systematic census for stellar transient events.

## Acknowledgments

We thank Yuichi Harikane and Anna Wright for helpful discussions. We thank Peter Williams for providing the three-color image of the HUDF used in Figure 4 and scripts for overplotting data. V.C. thanks the BlackHoleWeather project and PI Prof. Gaspari for salary support. V.C. also acknowledges support from the Chalmers Initiative on Cosmic Origins (CICO) as a visiting PhD student at Chalmers University and cofunding of his former PhD position at University of Trieste. J.C.T. acknowledges support from ERC advanced grant 788829 (MSTAR). A.R.Y. is supported by the Swedish National Space Agency (SNSA). M.J.H. is supported by the Swedish Research Council (Vetenskapsrådet) and is Fellow of the Knut & Alice Wallenberg Foundation. R.S.E. acknowledges generous financial support from the Peter and Patricia Gruber Foundation. M.J.H. and J.C.T. thank the staff in the Doors public house for service that germinated the ideas that led to this project.

*Facility: HST (WFC3).*

## Data Availability

All the HST data used in this paper can be found in MAST: DOI: 10.17909/7s5v-gz68.

# Appendix
## Variability in F105W and F160W

Figures 8 and 9 report the photometric variability as a function of the mean magnitude in the F105W and F160W filters, as detailed in Figure 1.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Photometric variability in the HUDF in the F105W filter from 2008–9 to 2012. Magnitude difference from epoch 1 (2008–9) to epoch 2 (2012) is plotted vs. average magnitude (gray points).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Photometric variability in the HUDF in the F160W filter from 2008–9 to 2012. Magnitude difference from epoch 1 (2008–9) to epoch 2 (2012) is plotted vs. average magnitude (gray points).</figcaption></figure>
ORCID iDs

Vieri Cammelli https://orcid.org/0000-0002-2070-9047

Jonathan C. Tan https://orcid.org/0000-0002-3389-9142

Alice R. Young https://orcid.org/0000-0001-9136-3701

Matthew J. Hayes https://orcid.org/0000-0001-8587-218X

Jasbir Singh https://orcid.org/0000-0002-6260-1165

Richard S. Ellis https://orcid.org/0000-0001-7782-7071

Aayush Saxena https://orcid.org/0000-0001-5333-9970

Nicolas Laporte https://orcid.org/0000-0001-7459-6335

Pierluigi Monaco https://orcid.org/0000-0003-2083-7564

Benjamin W. Keller https://orcid.org/0000-0002-9642-7193

References

Abel, T., Bryan, G. L., & Norman, M. L. 2002, *Sci*, 295, 93

Arnouts, S., Schiminovich, D., Ilbert, O., et al. 2005, *ApJL*, 619, L43

Bacon, R., Brinchmann, J., Conseil, S., et al. 2023, *A&A*, 670, A4

Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, *PASP*, 93, 5

Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, *Natur*, 553, 473

Banik, N., Tan, J. C., & Monaco, P. 2019, *MNRAS*, 483, 3592

Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, *AJ*, 132, 1729

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

Berg, D. A., Chisholm, J., Erb, D. K., et al. 2019, *ApJL*, 878, L3

Bertin, E., & Arnouts, S. 1996, *A&AS*, 117, 393

Bogdan, A., Goulding, A. D., Natarajan, P., et al. 2024, *NatAs*, 8, 126

Bouwens, R. J., Illingworth, G. D., Labbe, I., et al. 2011, *Natur*, 469, 504

Boyett, K., Bunker, A. J., Curtis-Lake, E., et al. 2019, *MNRAS*, 535, 1796

Bromm, V., Coppi, P. S., & Larson, R. B. 2002, *ApJ*, 564, 23

Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023, *A&A*, 677, A88

Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023, *A&A*, 677, A115

Cammelli, V., Monaco, P., Tan, J. C., et al. 2025, *MNRAS*, 536, 851

Casertano, S., de Mello, D., Dickinson, M., et al. 2000, *AJ*, 120, 2747

Chon, S., Hirano, S., Hosokawa, T., & Yoshida, N. 2016, *ApJ*, 832, 134

Cohen, S. H., Ryan, R. E. J., Straughn, A. N., et al. 2006, *ApJ*, 639, 731

den Brok, M., Seth, A. C., Barth, A. J., et al. 2015, *ApJ*, 809, 101

Ellis, R. S., McLure, R. J., Dunlop, J. S., et al. 2013, *ApJL*, 763, L7

Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, *AJ*, 132, 117

Finkelstein, S. L., & Bagley, M. B. 2022, *ApJ*, 938, 25

Freese, K., Gondolo, P., Sellwood, J. A., & Spolyar, D. 2009, *ApJ*, 693, 1563

Freese, K., Ilie, C., Spolyar, D., Valluri, M., & Bodenheimer, P. 2010, *ApJ*, 716, 1397

Freitag, M., Gürkan, M. A., & Rasio, F. A. 2006, *MNRAS*, 368, 141

Fruchter, A. S., & Hook, R. N. 2002, *PASP*, 114, 144

Graham, A. W. 2016, *ASSL*, 418, 263

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

Greif, T. H., & Bromm, V. 2006, *MNRAS*, 373, 128

Häberle, M., Neumayer, N., Seth, A., et al. 2024, *Natur*, 631, 285

Haehnelt, M. G., & Rees, M. J. 1993, *MNRAS*, 263, 168

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

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

Hayes, M. J., Tan, J. C., Ellis, R. S., et al. 2024, *ApJL*, 971, L16

Hirano, S., Hosokawa, T., Yoshida, N., et al. 2014, *ApJ*, 781, 60

Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011, *Sci*, 334, 1250

Hu, W., Papovich, C., Dickinson, M., et al. 2024, *ApJ*, 971, 21

Johnson, J. L., & Bromm, V. 2006, *MNRAS*, 366, 247

Kokorev, V., Caputi, K. I., Greene, J. E., et al. 2024, *ApJ*, 968, 38

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

Laporte, N., Nakajima, K., Ellis, R. S., et al. 2017, *ApJ*, 851, 40

Larson, R. L., Finkelstein, S. L., Kocevski, D. D., et al. 2023, *ApJL*, 953, L29

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

Lyu, J., Alberts, S., Rieke, G. H., & Rujopakarn, W. 2022, *ApJ*, 941, 191

Mainali, R., Kollmeier, J. A., Stark, D. P., et al. 2017, *ApJL*, 836, L14

Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2024, *A&A*, 691, A145

Matteo, T. D., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, *ApJ*, 676, 33

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

McKee, C. F., & Tan, J. C. 2008, *ApJ*, 681, 771

Mortlock, A., Conzelice, C. J., Hartley, W. G., et al. 2013, *MNRAS*, 433, 1185

Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, *Natur*, 474, 616

Mummery, A., & van Velzen, S. 2025, *MNRAS*, 541, 429

Natarajan, A., Tan, J. C., & O'Shea, B. W. 2009, *ApJ*, 692, 574

O'Brien, R., Jansen, R. A., Grogin, N. A., et al. 2024, *ApJS*, 272, 19

Pouliasis, E., Georgantopoulos, I., Bonanos, A. Z., et al. 2019, *MNRAS*, 487, 4285

Rafelski, M., Teplitz, H. I., Gardner, J. P., et al. 2015, *AJ*, 150, 31

Rieke, M. J., Robertson, B., Tacchella, S., et al. 2023, *ApJS*, 269, 16

Rindler-Daller, T., Montgomery, M. H., Freese, K., Winget, D. E., & Paxton, B. 2015, *ApJ*, 799, 210

Saxena, A., Pentericci, L., Mirabelli, M., et al. 2020, *A&A*, 636, A47

Schleicher, D. R. G., Reinoso, B., & Klessen, R. S. 2023, *MNRAS*, 521, 3972

Senchyna, P., Stark, D. P., Mirocha, J., et al. 2020, *MNRAS*, 494, 941

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

Singh, J., Monaco, P., & Tan, J. C. 2023, *MNRAS*, 525, 969

Spolyar, D., Freese, K., & Gondolo, P. 2008, *PhRvL*, 100, 051101

Stern, D., Eisenhardt, P., Gorjian, V., et al. 2005, *ApJ*, 631, 163

STSCI Development Team, 2012 DrizzlePac: HST Image Software, Astrophysics Source Code Library, record, ascl:1212.011

Suberlak, K. L., Ivezić, Ž., & MacLeod, C. 2021, *ApJ*, 907, 96

Sun, L., Wang, X., Teplitz, H. I., et al. 2024, *ApJ*, 962, 8

Susa, H., Hasegawa, K., & Tominaga, N. 2014, *ApJ*, 792, 32

Tan, J. C. 2008, in IAU Symp. 255, Low-Metallicity Star Formation: From the First Stars to Dwarf Galaxies, ed. L. K. Hunt, S. C. Madden, & R. Schneider (Cambridge: Cambridge Univ. Press), 24

Tan, J. C., Beltrán, M. T., Caselli, P., et al. 2014, *Protostars and Planets VI*, ed. H. Beuther, R.S. Klessen, C. P. Dullemond, & T. Henning (Tucson, AZ: Univ. Arizona Press), 149

Tan, J. C., & McKee, C. F. 2004, *ApJ*, 603, 383

Tan, J. C., Singh, J., Cammelli, V., et al. 2024, *arXiv*:2412.01828

Tan, J. C., Smith, B. D., & O'Shea, B. W. 2010, in AIP Conf. Ser. 1294, First Stars and Galaxies: Challenges for the Next Decade, 34, ed. D. J. Whalen, V. Bromm, & N. Yoshida (Melville, NY: AIP), 40

Tremmel, M., Karcher, M., Governato, F., et al. 2017, *MNRAS*, 470, 1121

Vika, M., Driver, S. P., Graham, A. W., & Liske, J. 2009, *MNRAS*, 400, 1451

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, *MNRAS*, 444, 1518

Wang, F., Yang, J., Fan, X., et al. 2021, *ApJL*, 907, L1

Wise, J. H., Regan, J. A., O'Shea, B. W., et al. 2019, *Natur*, 566, 85

Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, *ApJ*, 592, 645

Young, A. R., Hayes, M. J., Saldana-Lopez, A., et al. 2025, *arXiv*:2508.15905