# Glimmers in the Cosmic Dawn: A Census of the Youngest Supermassive Black Holes by Photometric Variability

Matthew J Hayes, Jonathan C Tan, Richard S Ellis, Alice R Young, Vieri Cammelli, Jasbir Singh, Axel Runnholm, Aayush Saxena, Ragnhild Lunnan, Benjamin W Keller, et al.

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

## To cite this version:

Matthew J Hayes, Jonathan C Tan, Richard S Ellis, Alice R Young, Vieri Cammelli, et al. Glimmers in the Cosmic Dawn: A Census of the Youngest Supermassive Black Holes by Photometric Variability. The Astrophysical Journal Letters, 2024, 971 (1), pp.L16. 10.3847/2041-8213/ad63a7. hal-04537846

HAL Id: hal-04537846

https://hal.science/hal-04537846v1

Submitted on 17 Jan 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>
<figure><img src="image_3.png" /></figure>

# Glimmers in the Cosmic Dawn: A Census of the Youngest Supermassive Black Holes by Photometric Variability*

Matthew J. Hayes¹, Jonathan C. Tan²,³, Richard S. Ellis⁴, Alice R. Young¹, Vieri Cammelli²,⁵,⁶, Jasbir Singh²,⁵,⁷, Axel Runnholm¹, Aayush Saxena⁸, Ragnhild Lunnan¹, Benjamin W. Keller⁹, Pierluigi Monaco⁵,⁶,¹⁰,¹¹, Nicolas Laporte¹², and Jens Melinder¹

¹ Stockholm University, Department of Astronomy and Oskar Klein Centre for Cosmoparticle Physics, AlbaNova University Centre, SE-10691 Stockholm, Sweden; matthew.hayes@astro.su.se

² Department of Space, Earth & Environment, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden

³ Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA

⁴ Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

⁵ Dipartimento di Fisica, Sezione di Astronomia, Università degli Studi di Trieste, via G.B. Tiepolo 11, I-34131 Trieste, Italy

⁶ INAF—Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, I-34131 Trieste, Italy

⁷ INAF—Osservatorio Astronomico di Brera, via Brera 28, I-20121 Milano, Italy

⁸ Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

⁹ Department of Physics and Materials Science, University of Memphis, 3720 Alumni Avenue, Memphis, TN 38152, USA

¹⁰ INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste TS, Italy

¹¹ IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy

¹² Aix Marseille Université, CNRS, CNES, LAM (Laboratoire d'Astrophysique de Marseille), UMR 7326, 13388 Marseille, France
*Received 2024 April 7; revised 2024 July 15; accepted 2024 July 16; published 2024 August 6*

## Abstract

We report the first results from a deep near-infrared campaign with the Hubble Space Telescope to obtain late-epoch images of the Hubble Ultra Deep Field, 10–15 yr after the first epoch data were obtained. The main objectives are to search for faint active galactic nuclei (AGN) at high redshifts by virtue of their photometric variability and measure (or constrain) the comoving number density of supermassive black holes (SMBHs), $n_{SMBH}$, at early times. In this Letter, we present an overview of the program and preliminary results concerning eight objects. Three variables are supernovae, two of which are apparently hostless with indeterminable redshifts, although one has previously been recorded as a $z \approx 6$ object precisely because of its transient nature. Two further objects are clear AGN at $z = 2.0$ and 3.2, based on morphology and/or infrared spectroscopy from JWST. Three variable targets are identified at $z = 6–7$ that are also likely AGN candidates. These sources provide a first measure of $n_{SMBH}$ in the reionization epoch by photometric variability, which places a firm lower limit of $3 \times 10^{-4} \text{ cMpc}^{-3}$. After accounting for variability and luminosity incompleteness, we estimate $n_{SMBH} \gtrsim 8 \times 10^{-3} \text{ cMpc}^{-3}$, which is the largest value so far reported at these redshifts. This SMBH abundance is also strikingly similar to estimates of $n_{SMBH}$ in the local Universe. We discuss how these results test various theories for SMBH formation.

*Unified Astronomy Thesaurus concepts: Active galactic nuclei (16); High-redshift galaxies (734)*

## 1. Introduction

The search for active galactic nuclei (AGN) in the early Universe proceeds for several reasons. Each counted AGN means at least one supermassive black hole (SMBH) is identified, which places a lower limit constraint on the abundance of SMBHs at early times. Accurate determination of this abundance may in turn be used to test SMBH formation theories (e.g., Rees 1978; Volonteri 2010; Banik et al. 2019; Inayoshi et al. 2020; Singh et al. 2023). The recovered number density and associated SMBH seeding mechanisms may in turn be used as an empirically motivated input to galaxy formation simulations, where AGN are fundamentally implicated in

galaxy evolution, luminosity functions, quenching, and other phenomena.

Some theoretical models (e.g., Banik et al. 2019) predict that large numbers of AGN exist in the reionization epoch, and recent observations from JWST may (at least in part) confirm the abundance to be larger than general community-wide expectations from before 2022 (e.g., Harikane et al. 2023; Larson et al. 2023; Maiolino et al. 2023). Most specifically, the families of “little red dots” at $z > 4$ (Furtak et al. 2023; Matthee et al. 2024; Kokorev et al. 2024; Wang et al. 2024) mark an unexpected discovery. These compact objects are identified at high-$z$ with a smoothly rising red spectral energy distribution (SED) and often showing broad Balmer emission lines. While it is not settled that the majority of these contain AGN, spectroscopic studies show that a substantial fraction may do so, which could imply large numbers of AGN in the reionization era. With AGN expected to exhibit high hydrogen ionizing photon production efficiencies and escape fractions, they could (once more) become implicated in reionization scenarios.

The challenge to this field comes from the fact that identifying AGN at the luminosity regimes of typical galaxies (below $L^*$, with $M_{UV} \gtrsim -20$) is observationally difficult. There

<figure><img src="image_4.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 526555. These observations are associated with programs 11563, 12498, and 17073.</footnote>

<footnote>Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
<table><caption>**Table 1**
Observing Epochs and Times</caption><thead><tr><th>Year</th><th>Filters</th><th>Orbits</th><th>GO No./PI</th></tr></thead><tbody><tr><td rowspan="2">2008–2009</td><td>F105W</td><td>24</td><td rowspan="2">11563/Illingworth</td></tr><tr><td>F160W</td><td>53</td></tr><tr><td rowspan="3">2012</td><td>F105W</td><td>72</td><td rowspan="3">12498/Ellis</td></tr><tr><td>F140W</td><td>30</td></tr><tr><td>F160W</td><td>26</td></tr><tr><td>2023</td><td>F140W</td><td>30</td><td>17073/Hayes</td></tr></tbody></table>

are numerous tracers of accretion-powered excitation in galaxies, including, but not limited to, X-ray or radio continuum emission that is too bright to be produced by star formation, the presence of broad emission lines, the presence of emission lines with ionization potentials too high to be explained by stellar photoionization or shocks, various ratios of emission line fluxes, and mid-infrared (MIR) colors consistent with dust heated to temperatures of ~1000 K. Yet the heterogeneous population of AGN typically satisfies subsets of these criteria, which become observationally more challenging/less complete at different rates as surveys target higher redshifts. This leads to SMBHs probably being undercounted, with potentially large numbers going unnoticed among the ostensibly star-forming galaxy population at high-z. Multiple diagnostics are mandatory to place stronger constraints on the AGN abundance.

Here we argue that the photometric variability that results from changes in the mass accretion rate of SMBHs can provide a completely independent and complementary probe of AGN. Monitoring for variability selects AGN from imaging data directly by phenomena related to the SMBH, without any biases of photometric preselection (color, luminosity, compactness, etc), and has been used successfully at intermediate redshifts in deep Hubble Space Telescope (HST) imaging (Cohen et al. 2006; Pouliasis et al. 2019; O’Brien et al. 2024). However, the highest redshifts require monitoring at infrared wavelengths in order to even detect galaxies: JWST would obviously be the best tool, but as such a recent observatory, it will take years to build the requisite time baselines when accounting for cosmological time dilation. HST provides a head start. WFC3 was mounted in 2009 and almost immediately began imaging some of the best-studied extra-galactic fields at $\lambda = 1-1.6$ µm. The deepest of these is the Hubble Ultra Deep Field (HUDF), which was imaged in 2009 and 2012 (Bouwens et al. 2010 and Ellis et al. 2013, respectively) and provides the optimal field for a search for faint AGN at high redshifts.

We reimaged the HUDF in 2023 to search for variable sources over a 10–15 yr time window. This is critical for very high-z systems where time dilation reduces this delay by a factor of $(1+z)$, to time spans of just ~1 yr in the rest frame. In this Letter, we present an overview of our survey and the observational data sets (Section 2) and some selected preliminary results (Section 3) that include supernovae, independently confirmed AGN at $z=2-3.5$, and newly discovered transients whose photometric/spectroscopic redshifts are consistent with sources in the reionization era. We provide some discussion regarding the SMBH abundance at $z>6$ in Section 4 and the implications for SMBH seeding mechanisms.

## 2. Observations, Data, and Methods

### 2.1. *Observations and Reductions*

In order to maximize our sensitivity to high-z variable AGN, we elected to target the HUDF because this field has provided samples at the highest comoving volume density. The field was first observed in the optical with the Advanced Camera for Surveys (ACS; Beckwith et al. 2006) and then with the critical near-infrared (NIR) mode with WFC3/IR in 2008–2009 under GO 11563 (UDF09; PI: Illingworth), which conducted imaging in three filters (F105W, F125W, and F160W) over 192 orbits. In 2012, the field was imaged again under GO 12498 (UDF12; PI: Ellis), which substantially deepened the F105W and F160W observations and added a fourth filter, F140W, to search for Lyman break galaxies at $z \sim 8$. With the singular goal of searching for photometric variability in all sources in the HUDF IR footprint, we reimaged the field in 2023 September in the F140W filter, matching the precise center, field orientation, and depth (30 orbits) of the GO 12498 observation.

We reduced the F140W image with the calwfc3 pipeline and astrodrizzle (STSCI Development Team 2012) software using the high-level science products (HLSP) obtained from the Mikulski Archive for Space Telescopes as reference images. We also rerReduced the F140W image from the UDF12 campaign to verify that our methods match the depth of the HLSP image. At the same time, we produced independent reductions of the F105W and F160W images at the 2009 and 2012 epochs to search for variable sources over the shorter, earlier time baseline. Full details of the observations and data reduction will be presented in V. Cammelli et al. (2024, in preparation) and A. Young et al. (2024, in preparation). This facilitates variability searches at three epochs: the span 2009–2012 is sampled by the F105W and F160W filters, while the duration 2012–2023 is sampled by the F140W filter only (see Table 1 for details).

### 2.2. *The Search for Variable Sources*

We searched for variables in each matched pair of filters using two techniques: comparison of the nuclear/central photometry of galaxies at different epochs and detection of residual sources in pair-subtracted images. We avoid regions toward the edge of the image where the dithering pattern has led to excess noise and consider only a central pointing of the HUDF covering $123'' \times 139''$. The photometric comparison follows a similar method to that of O’Brien et al. (2024): for each of the six WFC3/IR images, we identify sources using Source Extractor (Bertin & Arnouts 1996). In contrast to performing “global” galaxy photometry (e.g., with Kron-like or moment-centered apertures as Source Extractor would do), we seek to obtain photometry centered upon the brightest unresolved sources within each galaxy. We therefore rely only upon the brightest pixel coordinates reported by Source Extractor and the background images it produces as part of its background subtraction procedure. We first use astropy’s photutils package to produce precisely recentered coordinates in the immediate vicinity of the {XPEAK_IMAGE, YPEAK_IMAGE} for each source. At this coordinate, we then use photutils to perform aperture photometry in 4 pixel diameter (0.6"/26) apertures with a local background subtraction to remove local (noncompact) galaxy light. However, this would not include Poissonian uncertainties from the sky background, which we calculate from the background-level
images stored by Source Extractor (stored with CHECK-
IMAGE flags) using the same determined effective gains and
add in quadrature to the error per source using the local
photometry. We verified that this approach produces values
comparable to Source Extractor's MAG_APER measure-
ments but with more precise centroiding inside extended
sources.

For each of the six images we use for detection, we also
perform photometry in identical apertures in the other image
obtained in the same filter, similar to “dual image mode” in
Source Extractor. Thus, for each image, we obtain
photometry of each galaxy nucleus in detection and comparison
bands, working in both directions. We contrast these local
aperture magnitudes in each of the images with respect to the
photometric uncertainty of each source. However, we must also
account for the artificially underestimated uncertainties that
result from the signal being correlated between adjacent pixels
during the drizzling process (e.g., Casertano et al. 2000; Fruchter
& Hook 2002). We examine the global distribution on the Δm
versus m diagram (V. Cammelli et al. 2024, in preparation),
where Δm = m<sub>det</sub> − m<sub>test</sub>; here m<sub>det</sub> is the magnitude in the
image used for source detection, and m<sub>test</sub> is the magnitude in the
image used to test for photometric variability.

We first correct a ≈0.01 mag systematic offset that is
identified at all magnitudes, probably because of slightly
imperfect zero-points being reported to images obtained many
years apart. We then proceed by assuming that the over-
whelming majority of sources in the field are not variable. We
compute the one-sided standard deviation of the Δm distribu-
tion in bins of 0.5 mag, where we include only sources that are
brighter in m<sub>test</sub> compared to m<sub>det</sub> (hence the matched-aperture
photometry in both directions). We contrast the resulting
standard deviations with the propagated uncertainty on the
magnitude difference derived from our photometric catalogs:
σ<sub>m</sub> = (σ<sup>2</sup><sub>det</sub> + σ<sup>2</sup><sub>test</sub>)<sup>1/2</sup>. We then compute the total scale factor
that scales the combined uncertainty reported by astropy/
photutils to reproduce the global Δm distribution. This
ensures that, for example, 31.7% of our sources have
Δm/σ<sub>m</sub> ≥ 1, 4.5% of sources have Δm/σ<sub>m</sub> ≥ 2, etc. See
O'Brien et al. (2024) for a more detailed discussion.

To be identified as a variable source, an object must be
detected at Δm/σ<sub>m</sub> ≥ 2 in two or more filters or Δm/σ<sub>m</sub> ≥ 3 in
one filter. Should those two filters be the F105W and F160W
that were taken at the same epochs, we further demand that Δm
has the same sign (i.e., that the source is consistently inferred to
be brightening or fading over 2009–2012). If one of the filters
is F140W, we enforce no sign constraints on Δm, allowing for
a source to brighten and then fade, or vice versa. We have
verified that for the number of sources in the images, these
selections result in a negligible suspected level of contamina-
tion from nonvarying sources (V. Cammelli et al. 2024, in
preparation) at magnitudes as faint as 29. For the matched
sources, we compute their “total significance of variability” as
the quadrature addition of Δm/σ<sub>m</sub> in two bands (or three if the
source is variable in all filters).

In addition, we produced pair-subtracted images for each
filter to search for variable sources that may have been missed
by photometric comparison and to verify the compact nature of
variables. We first smoothed each image with a very small
Gaussian kernel of 0.5 pixels and directly subtracted images
obtained in the same filters. This produced six “difference
images.” We ran Source Extractor on the difference

images after recomputing the effective gains for each image, as
the subtraction increases the sky noise and mimics a shallower
integration time. Each source was inspected by eye, and
spurious artifacts were removed. We retain sources that have
5σ significance as variables.

## 2.3. Redshift Determination

We cross-correlate each of these sources against positions in
known redshift catalogs, including catalogs from the Very
Large Telescope (VLT)/MUSE (spectroscopic; Bacon et al.
2023), the JADES program using both NIRSpec for spectro-
scopic redshifts (Bunker et al. 2023) and NIRCam+HST for
photometric redshifts (Rieke et al. 2023), and the Ultraviolet
UDF (UVUDF) photometric catalogs (Rafelski et al. 2015). We
also cross-correlate with known AGN in the GOODS-S field,
compiled by Lyu et al. (2022). Further, we extract our own
photometry from all HST filters (11 bands from WFC3 and
ACS; Beckwith et al. 2006; Bouwens et al. 2010; Ellis et al.
2013; Koekemoer et al. 2013; Teplitz et al. 2013), JWST in
nine broadband filters from JADES (Rieke et al. 2023), and five
intermediate-band filters from the JEMS program (Williams
et al. 2023). Finally, we have reduced the deep slitless
spectroscopy obtained from JWST/NIRISS (Pirzkal et al.
2024) under the NGDEEP campaign (GO 2079; PI: Finkel-
stein) and extracted spectra at the location of all variable
sources, recentering on sources detected in NIRISS direct
imaging. See V. Cammelli et al. (2024, in preparation) for the
redshift distribution of the variable sources and A. Young et al.
(2024, in preparation) for an analysis of their morphologies and
photometric inferences.

3. Selected Variable Sources in the HUDF

In total, we can identify 71 sources with photometric
variability greater than 3σ over all epochs and filters. These
have redshifts in the range 0–7 but are concentrated at lower z.
We present preliminary results in this Letter, with the focus
upon a subset of interesting variable objects. These include two
confirmed AGN at z = 2–3.2, which demonstrate the variability
method in finding convincing AGN at redshifts where they can
be independently confirmed. The HUDF is arguably the best-
studied region of the extragalactic sky and has been targeted
across the electromagnetic spectrum: independent AGN
confirmation methods include deep optical and NIR spectrosc-
opy from VLT/MUSE and HST/grism, the deepest X-ray
imaging from Chandra, radio continuum imaging, and deep
Spitzer imaging to measure MIR colors (see Section 4.1 for
more details). We also show results for three transients that are
most likely supernovae, of which one in particular appears at
z ≈ 6 but is in fact an interloper at lower redshift. We finally
present our confident variable sources with redshifts in the
range z > 5–7.

3.1. *Variable AGN at Intermediate Redshifts*

To demonstrate the utility of photometric variability as an
AGN probe, we draw attention to objects 1051264 and
1052126 in Figure 1 and Table 2.

Object 1051264—This source lies in the center of an edge-
on disk galaxy. Image subtraction (Figure 1) shows the object
is pointlike. It decreased in brightness by 0.3 mag between
2009 and 2012 at very high significance (Table 2), although the
magnitude of this variation is small—(0.08 ± 0.025 mag) by
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Photometric variability of 1051264 at $z=2$ (upper panels) and 1052126 at $z=3.2$ (lower panels) sources. We show difference images of F105W and F140W to the left, where the earlier epoch is subtracted from the later. The color bar is symmetric around 0, so black objects get brighter with time and white objects become fainter. Panels to the right show the ACS/F606W image taken in 2005, averages of the F105W and F140W images over both epochs at which they were observed, and the JWST F277W image. Intensity images are calibrated in MJy $\text{Sr}^{-1}$ with color bars shown to the side. Cutouts are $2''4 \times 2''4$.</figcaption></figure>

<table><caption>**Table 2**

Properties of Two Example Intermediate-redshift AGN Detected by Variability</caption><thead><tr><th>ID</th><th>R.A.<br>J2000/Gaia</th><th>Decl.<br>J2000/Gaia</th><th>Mag 2009<br>F105W/F160W</th><th>Mag 2012<br>F105W/F160W</th><th>Mag 2012<br>F140W</th><th>Mag 2023<br>F140W</th><th>Redshift</th><th>Source</th></tr></thead><tbody><tr><td>1051264</td><td>53.180819</td><td>-27.787328</td><td>26.006 ± 0.013</td><td>26.332 ± 0.013</td><td>24.890 ± 0.015</td><td>24.909 ± 0.016</td><td>1.95</td><td>NIRISS</td></tr><tr><td>1052126</td><td>53.178480</td><td>-27.784036</td><td>24.493 ± 0.006</td><td>24.378 ± 0.003</td><td>24.092 ± 0.004</td><td>24.596 ± 0.006</td><td>3.190457 ± 0.000038</td><td>MUSE</td></tr></tbody></table>

global (Kron) photometry—because of dilution by the host galaxy. Its absolute magnitude is $M_{UV} = -20.3$, suggesting an AGN origin. Variability is confirmed in both F105W and F160W images, but no variability was detected in 2012–2023 in the F140W image. The source is obviously pointlike after the host is subtracted and is located at the nucleus of the galaxy, suggesting probable AGN activity.

The source is identified by Rafelski et al. (2015) and has a photometric redshift of 1.887. In Figure 2, we present the NIR spectra obtained using the F115W, F150W, and F200W filters of NIRISS. Emission lines of [O II], [O III], and Hα are unambiguous, which securely anchors the redshift to $z=1.95$, but there are no obvious signatures of AGN activity in this low-resolution spectrum.

**Object 1052126**—This object has a spectroscopic redshift of 3.1905 from VLT/MUSE. It is identified as a strong emitter of Lyα He II 1640, C IV 1550, and Si III 1883, as well as other high-ionization lines (Bacon et al. 2023), and is clearly an AGN. It is comparably bright, with $m_{F140} = 24.1$, and has varied by -0.53 and +0.50 mag from 2009 to 2012 and 2012 to 2023, respectively (brightening then dimming). The flux in the variable component alone amounts to an absolute magnitude of -21.3, corresponding approximately to the luminosity of an $L^*$ galaxy.

The complex of [O III] λλ4959, 5007 Å and Hβ is unambiguous in the NIRISS spectrum (Figure 2) at about 2.1 µm, anchoring the redshift to $z=3.2$, which is fully consistent with VLT/MUSE. Hγ emission is also clear in the F200W filter. The F150W filter shows an obvious line of [Ne III] λ3869, but more interestingly, there is a weak line at $\lambda_{\text{obs}} = 14400$ Å. This is identified as the [Ne V] λ3426 line, which requires a source of photons with $h\nu > 97$ eV or dense coronal gas with temperature $T > 10^5$ K. This can only be produced by an AGN. A strong Mg II feature is also identified in the F115W filter. Hβ is broader than the [O III] lines that are

<figure><img src="image_6.png" /><figcaption>**Figure 2.** JWST/NIRISS spectra of the sources shown in Figure 1. The three different wavelength regions show the F115W, F150W, and F200W filters, while the two spectra per wavelength range show extractions when the spectra are dispersed in the direction of rows or columns.</figcaption></figure>

observed in the same filter, with an FWHM of ≈4600 and ≈1400 km s⁻¹, respectively, where [O III] is consistent with instrumental broadening but Hβ is not. Hγ is of course weaker but appears to be as broad as Hβ, and Mg II shows a consistent FWHM. Most forbidden lines are too weak to provide reliable measurements of their velocity dispersion, but from the [O III] lines alone we see the familiar situation where permitted lines probe gas in the broad-line regions (BLRs) but BLR densities are sufficiently high that forbidden lines are collisionally
<figure><img src="image_7.png" /><figcaption>Figure 3. Same as Figure 1 for the probable supernovae.</figcaption></figure>

<table><caption>**Table 3**
Same as Table 2 but for the Supernova Candidate Sources</caption><thead><tr><th>ID</th><th>R.A.<br>J2000/Gaia</th><th>Decl.<br>J2000/Gaia</th><th>Mag 2009<br>F105W/F160W</th><th>Mag 2012<br>F105W/F160W</th><th>Mag 2012<br>F140W</th><th>Mag 2023<br>F140W</th><th>Redshift</th><th>Source</th></tr></thead><tbody><tr><td>1402129</td><td>53.154885</td><td>-27.793421</td><td>...</td><td>...</td><td>28.235 ± 0.042</td><td>26.575 ± 0.01</td><td>1.07<sup>+0.04</sup><sub>-0.02</sub></td><td>JADES phot</td></tr><tr><td>1402146</td><td>53.134812</td><td>-27.788965</td><td>...</td><td>...</td><td>≥30.0</td><td>27.234 ± 0.02</td><td>...</td><td>...</td></tr><tr><td>1051237</td><td>53.143873</td><td>-27.793327</td><td>26.927 ± 0.017</td><td>≥30.1</td><td>...</td><td>...</td><td>6.86</td><td>R15</td></tr></tbody></table>

deexcited, and [O III] emission probes only the more extended narrow line-emitting region.

## 3.2. Supernova Candidates

We show images of the three supernova candidates in Figure 3 and present general properties in Table 3.

**Object 1402129**—This is a clear stellar transient source exploding near the edge of a disk galaxy and is clearly not associated with nuclear activity. It is identified in our later epoch of F140W imaging and has increased in brightness by almost 2 mag—a 2012 magnitude is reported in Table 3, but this corresponds to the host galaxy in a matched aperture. Despite being relatively bright, the host galaxy only has a redshift reported from JADES photometry, which confirms the galaxy at $z \simeq 1.07$. With an absolute magnitude of $-16.9$ in the rest-frame R band, the supernova could have been of any type. Since this redshift would provide a factor of only $\simeq 2$ in cosmological time dilation, we expect the supernova to have been relatively recent: the slowly evolving superluminous supernovae (SLSNe) would typically have faded by 4 mag in 1 yr, with ordinary Type Ia and core-collapse supernovae fading faster still.

**1402146**—This is a remarkable transient source that was identified in our second epoch F140W imaging and exploded between 2012 and 2023. No host galaxy is apparent in the HST/IR imaging, but a faint source is just visible in the ACS imaging. Of our 15 visits in the observing program, one failed because of guide-star acquisition and was repeated 1 month

after the remainder of the observations. 1402146 was also visible in the repeat visit and had not moved, suggesting that it lies at cosmological distances and is not a solar system object. The source in the ACS image is not detected in redder bands and is so faint it does not enter any photometric redshift catalogs, so the distance to this source cannot be determined.

**1051237**—This object appears very similar to 1402146, except it had become fainter between the 2008 and 2012 images in the F105W and F160W filters, vanishing from 26.9 mag (both filters) to undetected. This source is intriguing, since it has a photometric redshift in the Rafelski et al. (2015) catalog of 6.84 and appears to be a supernova mistaken for a galaxy. 1051237 is visible in all the NIR filters of the HUDF09 campaign, and, while it is not visible in data from HUDF12, still appears in deep stacks that combine the two epochs. When the ACS optical data were obtained in 2005, the source had not yet exploded, but it had subsequently faded by the time the UVUDF (Teplitz et al. 2013) imaging was obtained. Consequently, in HST imaging only, the source is visible throughout the NIR but not in the UV/optical. In multiband photometry, it therefore appears as a z-band dropout and is assigned a best photometric redshift of $\sim 7$. The two galactic sources that lie about 1" to the north in the HST images have redshifts of 6, but other objects with redshifts in the range 0.2–1.1 lie within 2".

## 3.3. Three Variable Sources in the Reionization Era

After cross-correlating the coordinates of our variable objects with existing redshift catalogs (Section 2), we identify
<figure><img src="image_8.png" /><figcaption>Figure 4. Same as Figure 1 for the $z = 6-7$ sources.</figcaption></figure>

<table><caption>**Table 4**

General Properties of the $z > 6$ Variable Sources</caption><thead><tr><th>ID</th><th>R.A.<br>J2000/Gaia</th><th>Decl.<br>J2000/Gaia</th><th>Mag 2009<br>F105W/F160W</th><th>Mag 2012<br>F105W/F160W</th><th>Mag 2012<br>F140W</th><th>Mag 2023<br>F140W</th><th>Redshift</th><th>Source</th></tr></thead><tbody><tr><td>1052123</td><td>53.161915</td><td>-27.786990</td><td>28.308 ± 0.070</td><td>27.996 ± 0.028</td><td>27.595 ± 0.024</td><td>27.764 ± 0.029</td><td>6.74<sup>+0.073</sup><sub>-0.038</sub></td><td>JADES phot</td></tr><tr><td>1052156</td><td>53.154051</td><td>-27.766003</td><td>28.391 ± 0.067</td><td>28.181 ± 0.034</td><td>≥29.9</td><td>≥29.9</td><td>6.234</td><td>JADES spec</td></tr><tr><td>101159</td><td>53.160519</td><td>-27.785931</td><td>28.961 ± 0.366</td><td>≥30.1</td><td>≥29.9</td><td>≥29.9</td><td>6.54<sup>+2.45</sup><sub>-2.59</sub></td><td>JADES phot</td></tr></tbody></table>

three sources with significance greater than $3\sigma$ and best redshifts greater than 6. These are, naturally, much fainter than the sources identified at mid-$z$ above. They will be studied in detail in upcoming papers (V. Cammelli et al. 2024; A. Young et al. 2024; both in preparation), and here we present the three sources at $z > 6$ that are detected at better than $3\sigma$ in Figure 4 and Table 4. We also extracted NIRISS spectra from the NGDEEP survey but found no convincing detections of either continuum light or emission lines in any of the three. Given the faintness of these targets and the insensitivity of low-resolution spectroscopy to weak emission lines (expected in the UV), this is not surprising.

**Object 1052123**—This is our highest-significance variable source at $z > 6$. It is detected as variable in all three photometric bands and is visible in the pair-subtracted images of Figure 4. Its light curve is shown in Figure 5. It is pointlike and brightens between 2009 and 2012 (both F105W and F160W) by about 0.3 mag, with a very similar change in brightness detected at both wavelengths implying a constant color of the variable source. Interestingly, its F105W–F140W color is 0.42 mag, while the F140W–F160W color is 0.02 mag (both measured in 2012 epoch data). This is probably due to a Ly$\alpha$ break at 9400 Å, which produces a flux deficit in F105W leading to a photometric redshift of $z \approx 6.75$, while the continuum slope is flat at redder wavelengths ($\beta = -1.88 \pm 0.05$). The object then fades again to 2023, decreasing in brightness by 0.17 mag ($4.5\sigma$ significance).

The source is assigned a photometric redshift of 6.44 by Rafelski et al. (2015) and 6.74 in the JADES redshift catalog,

and we lean toward the higher value because the spectral break must lie within the F105W filter. Assuming this redshift and that the entire variable component can be attributed to an AGN, it would have an equivalent absolute magnitude of $-19.3$ AB. This is much less luminous than the variable AGN recovered at mid-$z$ (Section 3.1), which have $M_{UV} = -20.3$ and $-21.9$.

The whole duration of the light curve in Figure 5 is 15 yr, which corresponds to 2 yr in the rest frame, at filters probing $\lambda_{rest} \approx 1300-2000$ Å. These timescales exceed what is seen even for SLSNe, which are the slowest-evolving and most UV-bright supernova class known (Lunnan et al. 2018; Angus et al. 2019; Perley et al. 2020). In particular, while the peak UV luminosity is well within reach for an SLSN, the UV generally fades away on a timescale of approximately months due to cooling temperatures as the supernova expands (Yan et al. 2017). Similarly, while pair-instability supernovae are predicted to have very slow light curves due to their extremely high ejecta masses, they are also predicted to have red spectra thanks to the copious amounts of iron-peak elements produced in the explosion, which blankets the UV flux (Kasen et al. 2011; Dessart et al. 2013). We therefore favor an AGN over an extreme supernova interpretation.

**Object 1052156**—This object is identified as variable in the F105W and F160W filters (2009–2012) but is not found to be variable in the later epoch. The variability is clear in the pair-subtracted images, and if all the variability can be attributed to a single new source, its equivalent magnitude would be $-18.5$. This magnitude is also consistent with the object being a supernova, but given the change in brightness of 0.2 mag in
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Light curve for object 1052123. The lower x-axis shows the year of observation, while the upper x-axis shows the time since 2008 accounting for $(1+z)$. The left axis shows the observed apparent magnitude, while the right axis shows the absolute magnitude for $z=6.74$. Filters are shown in color.</figcaption></figure>

5 yr, it seems implausible (using the same argumentation as for 1052123). The combined significance across the two filters is only $3.5\sigma$ in this case, but even at $\approx 28$ mag there is still not an enormous parent sample of sources, and we do not expect to find one with a magnitude difference that is $3.5\sigma$ outside its distribution from pure chance. It is therefore probable that this source is a genuine variable. This source has a spectroscopic redshift of 6.23 (Bunker et al. 2023), as well as various photometric redshift estimates consistent with this.

**Object 101159**—The final object, 101159, is a remarkable variable source. It is clearly identified by the image subtraction method as a pointlike variable source that was visible in 2009 but not in 2012 or 2023. It is not visible in the ACS optical imaging obtained in 2006, implying a possible high redshift. It is clearly detected in the NIRCam/F277W image from JADES, which is the deepest in the JADES survey. The source likely comprises three diffuse clumps within the $1''$ region centered on the coordinates of the transient source, each of which is particularly faint, with magnitudes around $\sim 29$ in the F277W filter. The southern region is coincident with the location of the transient, which has no entry in the JADES catalogs. JADES reports a best-fit photometric redshift for the northwestern source of $z=6.54$, although the $1\sigma$ range is $3.84-8.99$. The northeastern source has a photometric redshift of $5.86$, consistent with the object nearest our transient source. If genuinely at $z>6$, the morphology of these sources is unlike most galaxies typically identified at high redshifts, where source identification is biased toward compact galaxies. This source has unfortunately decreased in brightness since 2008. The source could be an SLSN or an AGN, but in either case it outshines the host galaxy, which is only detectable by JWST at later epochs in redder images. We favor the AGN scenario for similar reasons to object 1052123: the probability of an SLSN explosion in our volume is $\sim 10^{-4}$, and these transients typically exhibit red colors and are hard to detect in the rest UV. We have verified that the source is visible in a deep stack of all the JADES imaging at NIR wavelengths and invisible in a stack of the ACS optical filters, again suggesting the object lies at high $z$.

**A note on interloping sources**—We initially identified another high-significance ($7\sigma$) variable at $\alpha, \delta = 53.1828095$,

$-27.7760521$, with $m_{F140W} = 28.1$. This object has a photometric redshift of 6.74 in the JADES catalog but is assigned $z_{\text{phot}} \simeq 1.7$ by Rafelski et al. (2015) and is clearly detected in ACS/F606W, implying that a low-$z$ solution must be correct.

# 4. Discussion

## 4.1. Variability as a Method to Find AGN and Estimate $n_{\text{SMBH}}$

We have identified several variable sources in the HUDF using photometric monitoring at three epochs: 2008, 2012, and 2023. In forthcoming papers (V. Cammelli et al. 2024, in preparation; A. Young et al. 2024, in preparation), we will present many more variable sources in the field at redshifts between 0 and 7. Here we have presented preliminary results for eight interesting targets, two of which are clear AGN at $z \simeq 2-3$, and three are probable supernovae that are also likely to have exploded at intermediate redshifts. We find three variable sources at redshifts greater than 6. Based upon their absolute magnitudes, which are measured at a rest wavelength of $\simeq 1300$ Å, we believe these variable sources to be AGN.

We argue that variability searches provide a valuable and efficient technique for finding AGN in deep imaging surveys. While only a fraction of AGN are detectable by variability, the advantage of monitoring is that it provides a complete, 100% multiplexed survey of everything in the field with no preselection. Furthermore, there is no additional challenge related to the different depths attainable for AGN diagnostics, as is the case for radio or X-ray surveys: if an object can be observed in imaging, it can be tested for variability.

Luminous quasars have been known at higher redshifts for some time (Fan et al. 2006; Mortlock et al. 2011; Bañados et al. 2018) and are all $\ge 4$ mag more luminous than our sources, and potentially higher-redshift AGN have also been reported (Harikane et al. 2023; Larson et al. 2023; Maiolino et al. 2023) at luminosities closer to $L^*$. However, as our survey examined only the volume of the HUDF and finds three probable AGN at $6<z<7$, this yields a higher comoving number density of AGN than previously reported.

Our survey samples a volume of $\simeq 10,500$ comoving Mpc$^3$ between redshifts 6 and 7, leading to a number density of $n_{\text{SMBH}} \sim 2.9 \times 10^{-4}$ cMpc$^{-3}$. However, only a fraction of AGN in a given luminosity range are expected to be detected by variability. Lyu et al. (2022) have recently published a sample of confirmed AGN in the GOODS-S field, within which the HUDF is entirely contained. Thirty-one of their AGN fall in the HUDF footprint that we target and have been identified by combinations of X-ray brightness, UV-to-MIR SED properties, optical spectral features, MIR colors, radio-loudness, continuum shape, and variability. We recover eight of these 31 AGN that have redshifts between 0.6 and 3.2 (object 1052162) by photometric variability. Lyu et al. (2022) identified these AGN by MIR colors (four), X-ray luminosities (seven), radio-loudness (seven), or optical spectroscopy (one). Some are obviously identified by more than one tracer. Importantly, none of our 71 prime candidates have been reported previously by photometric variability: this includes the HST imaging used in Lyu et al. (2022), who identified three variable AGN in the HUDF area, and Cohen et al. (2006), who reported 45 further sources. We note that all this imaging is shallower and obtained at shorter wavelengths and over shorter time baselines than our HUDF reimaging campaigns and should not necessarily match NIR long-baseline studies. None of our $z>6$ AGN are in the
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Cosmic evolution of SMBH abundance. The comoving number density of SMBHs, $n_{SMBH}$, is plotted vs. redshift, $z$. The observational constraint derived by counting $z=6–7$ AGN in the HUDF variability survey is shown by the red open square, and the lower bound shows the Poisson uncertainty. The red open diamond corrects for variability incompleteness, and the red filled diamond is for luminosity incompleteness (see text). Previous observational estimates at $z=0$ from local galaxies (Banik et al. 2019) and at $z=4–7$ from broad emission line sources (Harikane et al. 2023) are shown with the black square and diamonds, respectively. For the $z\sim4–7$ sources, we are unable to make any correction for the fraction of nonaccreting SMBHs and leave all points as lower limits. Various Pop III.1 formation models (Banik et al. 2019; 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) due to mergers. The green dashed line shows SMBH seeding based on an HMT above a mass of $7.1 \times 10^{10} M_{\odot}$ (Vogelsberger et al. 2014), with the shaded region illustrating a factor of 2 variation in this mass scale. The magenta dashed line shows results from a model of SMBH formation via “direct collapse” (Chon et al. 2016).</figcaption></figure>

Lyu et al. catalogs, almost certainly because they are too faint to be independently verified by the above techniques at current observational depths. See V. Cammelli et al. (2024, in preparation) for a more detailed discussion. From the ratio 8/31, we derive a crude estimate of this variability completeness correction factor of 3.85. After such a completeness correction is applied, $n_{SMBH} \sim 1.1 \times 10^{-3} \text{ cMpc}^{-3}$.

Next, it should be noted that our survey is only sensitive to AGN above a minimum luminosity of $M_{UV} \simeq -18.6$, so this estimate is still a lower limit on the true value of $n_{SMBH}$. Harikane et al. (2023) have detected AGN via the presence of broad emission lines and estimated $n_{SMBH}$ from $z \sim 4$ to 7 by extrapolating down an assumed luminosity function to $M_{UV} \simeq -17$. For a fair comparison of our results to those of Harikane et al. (2023), we carry out the same method and estimate a luminosity function (assuming $L^* = -20$ and $\alpha = -2$) correction factor of $\sim 7$, which boosts $n_{SMBH} \sim 7.7 \times 10^{-3} \text{ cMpc}^{-3}$. Note that to the extent that not all AGN are detected down to $M_{UV} \simeq -17$, this should still be regarded as a lower limit. In fact, at these redshifts, an empirical estimate of $n_{SMBH}$ can only be efficiently made by counting AGN and will include only SMBHs in a rapidly accreting state. The true comoving density would also need a final correction for the fraction of SMBHs that are currently dormant. The fact that we have no observational basis on which to estimate this final correction leads us to present all these densities as lower limits.

Our estimates of $n_{SMBH}$ at $z=6–7$ in the HUDF via variability and those of Harikane et al. (2023) via the presence of broad spectral lines from $z=4$ to 7 are shown in Figure 6. Here we also show a data point at $z=0$, i.e., $n_{SMBH}(z=0) \simeq 4.6 \times 10^{-3} \text{ cMpc}^{-3}$ with an uncertainty ranging from $10^{-3}$ to $10^{-2} \text{ cMpc}^{-3}$ (Banik et al. 2019). This is calculated by assuming

that each galaxy with luminosity greater than $L_{min} = 0.33L^*$ hosts an SMBH, with the error bar around this point assuming a range of $L_{min}$ from 0.1 to $1.0 L^*$ (see also Vika et al. 2009 for a similar estimated value of $8.79 \times 10^{-3} \text{ cMpc}^{-3}$). We see that our estimate of $n_{SMBH}$ in the HUDF is a factor of several higher than that of Harikane et al. (2023) at the equivalent redshifts but very similar to their estimate at $z \sim 4$ and the local estimate at $z=0$. Figure 6 also shows some predictions of theoretical models, which are discussed in the next section.

## 4.2. Implications for SMBH Formation Theories

Our estimate of $n_{SMBH}(z=6–7) \gtrsim 8 \times 10^{-3} \text{ cMpc}^{-3}$ has implications for SMBH seeding schemes and formation theories. Many implementations of SMBH seeding in cosmological simulations use simple threshold conditions. For example, a “halo mass threshold” (HMT) model was used in the Illustris and Illustris-TNG simulations (e.g., Vogelsberger et al. 2014), which seeded SMBHs in halos reaching a mass of $7.1 \times 10^{10} M_{\odot}$. The abundance of such massive halos and associated SMBHs is shown as a function of redshift in Figure 6. By $z \simeq 6.5$, this model has $n_{SMBH} \sim 10^{-3} \text{ cMpc}^{-3}$ but with its value rising rapidly as the Universe evolves to lower redshifts. By $z=0$, it has reached $n_{SMBH} \simeq 2 \times 10^{-2} \text{ cMpc}^{-3}$, which is marginally consistent with observational estimates $n_{SMBH}(z=0) \sim 10^{-3}–10^{-2} \text{ cMpc}^{-3}$ (Vika et al. 2009; Banik et al. 2019). However, note that if there is a significant population of SMBHs that have been ejected from their host galaxies and/or have very low luminosities, then $n_{SMBH}(z=0)$ may have been underestimated.

“Direct collapse” is a physical model of SMBH formation from metal-free, UV-irradiated, atomically cooled (i.e., relatively massive) halos that has been proposed as a way to form
massive seeds early in the Universe (e.g., Bromm & Loeb 2003). However, simulations of this mechanism have struggled to produce sufficient numbers of SMBHs compared to the known $z=0$ population. For example, Chon et al. (2016) found $n_{\text{SMBH}} \sim 10^{-4} \text{cMpc}^{-3}$, while Wise et al. (2019) estimated a global $n_{\text{SMBH}} \sim 10^{-7}-10^{-6} \text{cMpc}^{-3}$.

The “Pop III.1” SMBH formation model (Banik et al. 2019; Singh et al. 2023) invokes the physical mechanism of dark matter annihilation changing the structure of primordial protostars (Spolyar et al. 2008; Natarajan et al. 2009; Rindler-Daller et al. 2015) allowing efficient accretion of the baryonic content of their minihalos to become supermassive stars with masses of $\sim 10^5 M_{\odot}$, which then collapse to SMBHs. However, only the first minihalos to form in each local region of the Universe, i.e., being pristine and undisturbed by external feedback, are Pop III.1 sources and undergo this evolution. The vast majority of minihalos are Pop III.2, i.e., irradiated by UV and/or shocked by supernova blast waves leading to enhanced hydrogen-deuterium cooling and formation of relatively low-mass stars with $\sim 10 M_{\odot}$ (Greif & Bromm 2006; Johnson & Bromm 2006). Banik et al. (2019) and Singh et al. (2023) have predicted the evolution of the comoving number density of SMBHs forming from Pop III.1 seeds. The main parameter of the Pop III.1 model is the “isolation distance” ($d_{\text{iso}}$), defined as the minimum separation distance a Pop III.1 source needs to have from already formed sources, with this expected to be set by the physical mechanism of radiative feedback. To reproduce the known population of low-$z$ SMBHs with $n_{\text{SMBH}}(z=0) \sim 5 \times 10^{-3} \text{cMpc}^{-3}$ requires $d_{\text{iso}} \simeq 100 \text{ kpc}$ (proper distance), i.e., $\sim 3 \text{ cMpc}$ at the typical epoch of formation at $z \sim 30$. These models predict that $n_{\text{SMBH}}$ is approximately constant with redshift from $z \sim 25$ down to $z=0$, with only modest decreases due to mergers (see Figure 6).

The HUDF variability estimate of $n_{\text{SMBH}}(z=6-7) \gtrsim 8 \times 10^{-3} \text{cMpc}^{-3}$ places constraints on $d_{\text{iso}}$ to be $\lesssim 100 \text{ kpc}$. For HMT models, it requires the HMT for SMBH formation to be $\lesssim 3 \times 10^{10} M_{\odot}$. However, such models begin to have more severe tension by having significantly greater abundance than the $z=0$ estimate of $n_{\text{SMBH}}$. Finally, the HUDF variability estimate of $n_{\text{SMBH}}$ is about $100 \times$ greater than the direct-collapse prediction of Chon et al. (2016) and at least $10^4 \times$ greater than that of Wise et al. (2019).

Another test of the Pop III.1 model that can be made from our detected AGN is to examine the minimum separation distance between AGN. Note that the size of the field of view of the HUDF at $z \sim 6$ is about 5.5 cMpc. The closest pair of high-$z$ AGN is sources 101159 and 1052123. The former has a larger uncertainty in its redshift. If these sources are at the same redshift of $z \simeq 6.7$, then their plane-of-sky separation corresponds to about 0.25 cMpc. Thus, we see that the Pop III.1 model with $d_{\text{iso}} \simeq 100 \text{ kpc}$ (proper), corresponding to a comoving separation of about 3 cMpc (and for negligible relative motion by $z \sim 7$), implies that there must be a redshift difference of $\Delta z \gtrsim 0.01$ between these AGN. Future spectroscopic measurements of the redshifts of these sources can test this prediction. Similarly, another prediction of the Pop III.1 model is that all these high-$z$ AGN are powered by single SMBHs, with binary AGN only emerging at lower redshifts (J. Singh et al. 2024, in preparation).

### 4.3. Future Prospects for Variability Searches

It is clear from the above analysis that larger samples of AGN at high $z$ are required to reduce the statistical (i.e.,

Poisson) and systematic uncertainties in the estimates of $n_{\text{SMBH}}$. JWST is required to push to detection of fainter AGN via variability, requiring an investment similar to that outlined in Jha et al. (2024). However, it will take years to establish monitoring campaigns that span significant time baselines, especially when accounting for the factor $(1+z)$ in cosmological time dilation. In contrast, HST’s legacy of deep NIR imaging already stretches back $\simeq 15$ yr, providing an excellent baseline for monitoring. The deeper IR regions of the CANDELS fields (Grogin et al. 2011; Koekemoer et al. 2011) of strong gravitational lensing clusters Lotz et al. (2017) would be the most appropriate next targets to carry out high-$z$ AGN variability searches to build more stringent constraints on theoretical models of SMBH formation.

## Acknowledgments

We thank Yuichi Harikane, Marta Volonteri, and Naoki Yoshida for helpful discussions and the constructive and thoughtful report from the anonymous referee that has improved the content of the manuscript. M.H. acknowledges the support of the Swedish Research Council, Vetenskapsrådet, and is Fellow of the Knut and Alice Wallenberg Foundation. J.C.T. acknowledges support from ERC Advanced grant 788829 (MSTAR). A.Y. is funded by the Swedish National Space Agency, SNSA. R.L. is supported by the European Research Council (ERC) under the European Union’s Horizon Europe research and innovation program (grant agreement No. 10104229—Tran-sPIre). We thank the staff in The Doors public house for service that germinated the ideas that led to this project. All the HST data used in this Letter can be found in MAST:10.17909/7s5v-gz68.

Facilities: HST (WFC3), JWST (NIRISS).

## ORCID iDs

Matthew J. Hayes https://orcid.org/0000-0001-8587-218X
Jonathan C. Tan https://orcid.org/0000-0002-3389-9142
Richard S. Ellis https://orcid.org/0000-0001-7782-7071
Alice R. Young https://orcid.org/0000-0001-9136-3701
Vieri Cammelli https://orcid.org/0000-0002-2070-9047
Jasbir Singh https://orcid.org/0000-0002-6260-1165
Axel Runnholm https://orcid.org/0000-0002-1025-7569
Aayush Saxena https://orcid.org/0000-0001-5333-9970
Ragnhild Lunnan https://orcid.org/0000-0001-9454-4639
Benjamin W. Keller https://orcid.org/0000-0002-9642-7193
Pierluigi Monaco https://orcid.org/0000-0003-2083-7564
Nicolas Laporte https://orcid.org/0000-0001-7459-6335
Jens Melinder https://orcid.org/0000-0003-0470-8754

## References

Angus, C. R., Smith, M., Sullivan, M., et al. 2019, MNRAS, 487, 2215

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

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

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

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2010, ApJL, 709, L133

Bromm, V., & Loeb, A. 2003, ApJ, 596, 34

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

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
Dessart, L., Waldman, R., Livne, E., Hillier, D. J., & Blondin, S. 2013, MNRAS, 428, 3227

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

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

Furtak, L. J., Zitrin, A., Plat, A., et al. 2023, ApJ, 952, 142

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

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

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

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

Jha, S. W., Casetti-Dinescu, D. I., Bernstein, G. M., et al. 2024, arXiv:2405.12297

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

Kasen, D., Woosley, S. E., & Heger, A. 2011, ApJ, 734, 102

Koekemoer, A. M., Ellis, R. S., McLure, R. J., et al. 2013, ApJS, 209, 3

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

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

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

Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97

Lunnan, R., Chornock, R., Berger, E., et al. 2018, ApJ, 852, 81

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

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

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

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

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

Perley, D. A., Fremling, C., Sollerman, J., et al. 2020, ApJ, 904, 35

Pirzkal, N., Rothberg, B., Papovich, C., et al. 2024, ApJ, 969, 90

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

Rees, M. J. 1978, Obs, 98, 210

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

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

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

STSCI Development Team, 2012 Astrophysics Source Code Library, ascl:1212.011

Teplitz, H. I., Rafelski, M., Kurczynski, P., et al. 2013, AJ, 146, 159

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

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

Wang, B., de Graaff, A., Davies, R. L., et al. 2024, arXiv:2403.02304

Williams, C. C., Tacchella, S., Maseda, M. V., et al. 2023, ApJS, 268, 64

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

Yan, L., Quimby, R., Gal-Yam, A., et al. 2017, ApJ, 840, 57