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

# GLIMPSE: An Ultrafaint $\simeq 10^5 M_\odot$ Pop III Galaxy Candidate and First Constraints on the Pop III UV Luminosity Function at $z \simeq 6-7$

Seiji Fujimoto<sup>1,2,3,28</sup> (D), Rohan P. Naidu<sup>4,27</sup> (D), John Chisholm<sup>3</sup> (D), Hakim Atek<sup>5</sup> (D), Ryan Endsley<sup>3</sup> (D), Vasily Kokorev<sup>3</sup> (D), Lukas J. Furtak<sup>6</sup> (D), Richard Pan<sup>7</sup> (D), Boyuan Liu<sup>8</sup> (D), Volker Bromm<sup>3</sup> (D), Alessandra Venditti<sup>3</sup> (D), Eli Visbal<sup>9</sup> (D), Richard Sarmento<sup>10</sup> (D), Andrea Weibel<sup>11</sup> (D), Pascal A. Oesch<sup>11,12,13</sup> (D), Gabriel Brammer<sup>14</sup> (D), Daniel Schaerer<sup>11</sup> (D), Angela Adamo<sup>15</sup> (D), Danielle A. Berg<sup>3</sup> (D), Rachel Bezanson<sup>16</sup> (D), Rychard Bouwens<sup>17</sup> (D), Iryna Chemerynska<sup>5</sup> (D), Adélaïde Claeyssens<sup>15</sup> (D), Miroslava Dessauges-Zavadsky<sup>11</sup> (D), Anna Frebel<sup>4</sup> (D), Damien Korber<sup>11</sup> (D), Ivo Labbe<sup>18</sup> (D), Rui Marques-Chaves<sup>11</sup> (D), Jorryt Matthee<sup>19</sup> (D), Kristen B. W. McQuinn<sup>20,21</sup> (D), Julian B. Muñoz<sup>3</sup> (D), Priyamvada Natarajan<sup>22,23,24</sup> (D), Alberto Saldana-Lopez<sup>25</sup> (D), Katherine A. Suess<sup>26</sup> (D), Marta Volonteri<sup>5</sup> (D), and Adi Zitrin<sup>6</sup> (D)

<sup>1</sup>David A. Dunlap Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada; seiji.fujimoto@utoronto.ca

<sup>2</sup>Dunlap Institute for Astronomy and Astrophysics, 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada

<sup>3</sup>Department of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA

<sup>4</sup>MIT Kavli Institute for Astrophysics and Space Research, 70 Vassar Street, Cambridge, MA 02139, USA; rnaidu@mit.edu

<sup>5</sup>Institut d'Astrophysique de Paris, CNRS, Sorbonne Université, 98bis Boulevard Arago, 75014, Paris, France

<sup>6</sup>Department of Physics, Ben-Gurion University of the Negev, P.O. Box 653, Be'er-Sheva 84105, Israel

<sup>7</sup>Department of Physics & Astronomy, Tufts University, MA 02155, USA

<sup>8</sup>Institut für Theoretische Astrophysik, Zentrum für Astronomie, Universität Heidelberg, D-69120 Heidelberg, Germany

<sup>9</sup>Department of Physics and Astronomy and Ritter Astrophysical Research Center, University of Toledo, 2801 W Bancroft Street, Toledo, OH 43606, USA

<sup>10</sup>School of Earth and Space Exploration, Arizona State University, P.O. Box 871404, Tempe, AZ, 85287-1404, USA

<sup>11</sup>Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland

<sup>12</sup>Cosmic Dawn Center (DAWN), Copenhagen, Denmark

<sup>13</sup>Niels Bohr Institute, University of Copenhagen, Jagtvej 128, København N, DK-2200, Denmark

<sup>14</sup>Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen, Jagtvej 128, København N, DK-2200, Denmark

<sup>15</sup>Department of Astronomy, Oskar Klein center, Stockholm University, AlbaNova University center, SE-106 91 Stockholm, Sweden

<sup>16</sup>Department of Physics and Astronomy and PITT PACC, University of Pittsburgh, Pittsburgh, PA 15260, USA

<sup>17</sup>Leiden Observatory, Leiden University, NL-2300 RA Leiden, The Netherlands

<sup>18</sup>Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Melbourne, VIC 3122, Australia

<sup>19</sup>Institute of Science and Technology Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria

<sup>20</sup>Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

<sup>21</sup>Department of Physics & Astronomy, Rutgers, The State University of New Jersey, Piscataway, NJ 08854, USA

<sup>22</sup>Department of Astronomy, Yale University, 219 Prospect Street, New Haven, CT 06511, USA

<sup>23</sup>Department of Physics, Yale University, 217 Prospect Street, New Haven, CT 06511, USA

<sup>24</sup>Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA

<sup>25</sup>Department of Astronomy, Oskar Klein Centre, Stockholm University, 106 91 Stockholm, Sweden

<sup>26</sup>Department for Astrophysical & Planetary Science, University of Colorado, Boulder, CO 80309, USA

Received 2025 January 21; revised 2025 June 11; accepted 2025 June 16; published 2025 August 4

## Abstract

Detecting the first generation of stars, Population III (Pop III), has been a long-standing goal in astrophysics, yet they remain elusive even in the JWST era. Here we present a novel NIRCam-based selection method for Pop III galaxies, and carefully validate it through completeness and contamination simulations. We systematically search $\simeq 500$ arcmin<sup>2</sup> across JWST legacy fields for Pop III candidates, including GLIMPSE, which, assisted by gravitational lensing, has produced JWST's deepest NIRCam imaging thus far. We discover one promising Pop III galaxy candidate (GLIMPSE-16043) at $z = 6.50_{-0.24}^{+0.03}$, a moderately lensed galaxy ($\mu = 2.9_{-0.2}^{+0.1}$) with an intrinsic UV magnitude of $M_{UV} = -15.89_{-0.14}^{+0.12}$. It exhibits key Pop III features: strong H$\alpha$ emission (rest-frame EW 2810 $\pm$ 550 Å); a Balmer jump; no dust (UV slope $\beta = -2.34 \pm 0.36$); and undetectable metal lines (e.g., [O III]; [O III]/H$\beta < 0.44$), implying a gas-phase metallicity of $Z_{gas}/Z_\odot < 0.5\%$. These properties indicate the presence of a nascent, metal-deficient young stellar population ($<5$ Myr) with a stellar mass of $\simeq 10^5 M_\odot$. Intriguingly, this source deviates significantly from the extrapolated UV-metallicity relation derived from recent JWST observations at $z = 4-10$, consistent with UV enhancement by a top-heavy Pop III initial mass function or the presence of an extremely metal-poor active galactic nucleus. We also derive the first observational constraints on the Pop III UV luminosity function at $z \simeq 6-7$. The volume density of GLIMPSE-16043 ($\simeq 10^{-4} \text{ cMpc}^{-3}$) is in excellent agreement with theoretical predictions, independently reinforcing its plausibility. This study demonstrates the power of our novel NIRCam method to finally reveal

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

<footnote n="27"><sup>27</sup> These authors contributed equally to this work.</footnote>

<footnote n="28"><sup>28</sup> Hubble Fellow.</footnote>
distant galaxies even more pristine than the Milky Way's most metal-poor satellites, thereby promising to bring us closer to the first generation of stars than we have ever been before.

*Unified Astronomy Thesaurus concepts: Population III stars (1285); Early universe (435); High-redshift galaxies (734); Young star clusters (1833); Interstellar medium (847); Photoionization (2060); Reionization (1383); Emission nebulae (461); Galaxy formation (595); Galaxy evolution (594); Intermediate-mass black holes (816)*

# 1. Introduction

Members of the first generation of stars that formed from primordial, i.e., metal-free, gas are referred to as Population III (Pop III) stars (e.g., V. Bromm 2013; A. Frebel & J. E. Norris 2015; K. Inayoshi et al. 2020; R. S. Klessen & S. C. O. Glover 2023). Hypothesized since the 1960s (e.g., T. Page 1966), detecting any signatures of the presumably very massive and short-lived Pop III stars remains one of the elusive, ultimate frontiers of observational astrophysics. In many ways, Pop III stars represent the beginning of our beginnings. By forging metals for the very first time, these stars enriched the sterile ocean of hydrogen and helium gas left behind by the Big Bang, setting in motion the chemical evolution of the Universe that would one day lead to life on Earth. If not for their legacy of metals, the hot gas entrained in dark matter halos would struggle to cool and condense into stars and galaxies. And if not for the compact remnants they left behind, the behemoth supermassive black holes rippling spacetime today may have never been seeded. Few goals in astrophysics seem more fundamental than discovering these stars that are the very first chapter of our origin story.

Detecting light from the first Pop III stars forming at $z \approx 15-20$ demands imaging depths of $\approx 40$ mag, which is beyond the reach of all current and planned future facilities (e.g., A. T. P. Schauer et al. 2020). However, recent simulations suggest that Pop III star formation may have persisted for some time, down to lower redshifts, $z \simeq 6-7$ (e.g., C. Scannapieco et al. 2005; R. Jimenez & Z. Haiman 2006; L. Tornatore et al. 2007; A. Pallottini et al. 2014a; H. Katz et al. 2023; A. Venditti et al. 2023, 2024b). In these models, pristine pockets of gas that give rise to Pop III star formation in small dark matter host halos survive unpolluted due to, e.g., their remote locations on the edges of the then-emerging cosmic web and forming structures, inefficient feedback, and metal mixing in their immediate surroundings, or their residence in low-mass dark matter halos that have not yet undergone their first episode of star formation. When considering an observational search for early galactic systems at high redshift that contain Pop III stars, low-mass, faint ($L < 1\%L^*$) galaxies are hence the prime targets. As we demonstrate in Section 2.6, a pure Pop III-dominated source forming stars from gas in a primordial, metal-free dark matter halo is expected to have a stellar mass of $\approx 10^5 M_\odot$ requiring a depth of $\approx 31$ mag for detection.

What might the smoking gun signatures of Pop III stars in such systems be? Made from just H and He, massive Pop III stars embedded in primordial H II regions should produce spectra revealing strong hydrogen (e.g., Hα, Hβ) and helium lines (e.g., He II λ1640, 4686; J. Tumlinson et al. 2001; D. Schaerer 2002) accompanied by a complete absence of metal lines (e.g., [O III] λ5008). For pristine gas or extremely metal-poor conditions, top-heavy initial mass functions (IMFs) are expected, perhaps including extremely massive stars reaching $\approx 100 M_\odot - 1000 M_\odot$ (e.g., S. Hirano et al. 2014; S. Chon & K. Omukai 2025). The consequent hot and highly

ionized ISM conditions of primordial gas are expected to result in a steep Balmer jump in the nebular continuum portion of the spectrum (e.g., E. Zackrisson et al. 2011).

Prior to JWST (J. P. Gardner et al. 2023), efforts to search for Pop III stars at $z \gtrsim 6$ were limited to identifying rest-frame-UV signatures (Lyα and He II; e.g., D. Sobral et al. 2015), typically in fairly luminous systems ($\gtrsim 0.1L^*$) in order to be accessible to ground-based spectroscopy. JWST's sensitivity now allows us to push into the crucial low-mass dwarf galaxy regime where Pop III-dominated systems are most likely to occur, while also providing access to the full suite of rest-frame-UV to rest-frame-optical diagnostics. The question is therefore, what is the most efficient way of deploying JWST to scan a large number of sources for rare Pop III signatures?

Recent NIRSpec IFU results are tantalizing, identifying a candidate He II emitting clump at $z = 10.6$ (R. Maiolino et al. 2024b) and “LAP1,” a highly magnified [O III]-faint source at $z = 6.6$ (E. Vanzella et al. 2023a). Intriguingly, these source detections lack stellar continuum in their spectra, raising the possibility that instead of metal-poor stars, these two systems may instead represent nebular emission from pristine gas clouds due to irradiation by nearby star-forming regions (or an active galactic nucleus, AGN). However, given the limited survey volume (field of view of $3'' \times 3''$), the IFU, and more broadly, slit-spectroscopy, offer a more viable approach for follow-up studies (e.g., A. Venditti et al. 2024b).

NIRCam grism spectroscopy captures every source in the field of view, offering a combination of superb statistics, high-resolution (e.g., to split [O III]+Hβ out to $z \approx 9$), and spatially resolved spectra to identify metal-poor clumps (e.g., D. Kashino et al. 2023; J. Matthee et al. 2024; P. A. Oesch et al. 2023). However, this comes at the cost of sensitivity, as a result of which the deepest grism surveys are yet to push below $\approx 29$ mag. This implies that meaningful Pop III searches with the grism hinge on the discovery of extremely magnified arcs (e.g., R. P. Naidu et al. 2024).

The situation at hand is rather surprising—even after a full 2 yr of JWST operations, only a single source (LAP1) has been robustly confirmed below the “metallicity floor” of the Milky Way’s ultrafaint dwarf galaxy satellite population ([Fe/H]~ −2.5; e.g., S. W. Fu et al. 2023). Intriguingly, these local satellite galaxies may bear signatures of Pop III stars that exploded at $z > 6$ (e.g., M. Jeon et al. 2017). Remarkably, then, the best examples of primordial stellar populations are still at $z \sim 0$ and not at $z > 6$. This is despite the fact that these surviving systems formed their near-pristine metal-poor stars at exactly the redshifts ($z \gtrsim 6$) that are now well within our grasp with JWST (e.g., T. M. Brown et al. 2014; D. R. Weisz et al. 2014; A. Savino et al. 2023). At these early epochs, however, JWST is yet to push to the low-metallicity regimes and low stellar masses that the local dwarf satellite galaxies span. Ultradeep observations and new selection techniques are required to provide a breakthrough if we ever want to capture
the high-redshift ancestors of today's ancient ultrafaint satellite galaxies.

Indeed, an efficient approach that is yet to be systematically exploited involves capturing Pop III signatures through deep imaging surveys that deploy strategically chosen medium bands and broad bands. While the maximal rest-frame equivalent width (EW) of the rest-UV helium lines is $\simeq 100-200$ Å (e.g., D. Schaerer 2002; K. Nakajima & R. Maiolino 2022), the strong rest-frame optical $H\alpha$ line reaches EWs of $\gtrsim 2000-3000$ Å in young stellar populations (e.g., R. Endsley et al. 2023). This could significantly boost the flux even in broadband filters. In fact, this flux excess in broadband photometry has been used to constrain the EWs of the strong rest-frame optical $H\alpha$ and [O III] lines of $z \simeq 6-9$ galaxies with the Hubble Space Telescope (HST) and Spitzer for several years (e.g., I. Labbé et al. 2013; R. Smit et al. 2015; G. W. Roberts-Borsani et al. 2016; M. Stefanon et al. 2022). The same experiment may now be performed with JWST/NIRCam out to $z \simeq 9$ (e.g., R. Endsley et al. 2023; M. Llerena et al. 2024; C. Simmonds et al. 2024) and even at $z > 10$ with MIRI, albeit at much reduced efficiency (e.g., J. M. Helton et al. 2025). This provides a new opportunity for finding Pop III star systems, provided the classical line-flux excess analysis is optimized to capture the unique SED shapes of Pop III-dominated galaxies exhibiting strong $H\alpha$ and $H\beta$ but lacking any [O III] lines.

The redshift range of $z \simeq 6-7$ presents an ideal window to identify Pop III galaxies purely with NIRCam using the extensive rest-frame UV and optical information available. This is the highest redshift where $H\alpha$ is accessible. Furthermore, NIRCam's sensitivity peaks at $3.5-4$ µm (J. R. Rigby et al. 2025). This implies that the most stringent upper limit on the [O III] line at 5008 Å may be placed at $z \sim 6.5$. Higher redshifts may be probed using NIRCam +MIRI (E. Zackrisson et al. 2011; J. A. A. Trussler et al. 2023; see also T. Nagao et al. 2008; A. K. Inoue 2011), but MIRI's smaller field of view and much lower throughput provide significant overall limitations. Given the large volumes already surveyed and planned with NIRCam, a pure NIRCam search at $z \simeq 6-7$ not only lays a promising foundation for the discovery of early Pop III galaxies, but also enables the provision of stringent constraints on their volume density. Importantly, even nondetections provide valuable insights, offering strong constraints on Pop III star formation predictions, enrichment processes in chemically unevolved dark matter halos, and feedback mechanisms in early galaxies. For example, even upper limits on the Pop III star formation rate at these redshifts would guide the Roman mission in undertaking dedicated surveys for the putative pair-instability supernovae (PISNe), whose yields remain speculative at present (e.g., T. J. Moriya et al. 2022). Similarly, forecasts for the number of events (e.g., black hole mergers) detected by future gravitational-wave observatories (e.g., LISA) critically depend on the still observationally unconstrained Pop III IMF.

In this paper, we present our novel technique to select the best candidate Pop III galaxies using JWST/NIRCam data alone, and the first systematic search of Pop III galaxies in publicly available data from deep JWST legacy surveys. Notably, our search includes the Cycle 2 large program GLIMPSE (PID 3293; PIs H. Atek & J. Chisholm) that has produced JWST's deepest NIRCam imaging data set to date (accounting for lensing). It allows, for the first time, to push

into the early low-mass dwarf galaxy regime where Pop III-dominated galaxies are the likeliest to occur. By leveraging strong lensing produced by the well-studied massive galaxy cluster of Abell S1063 (D. Lutz et al. 2016), GLIMPSE is designed to detect ultrafaint galaxies in the epoch of reionization. The GLIMPSE filter set is composed of seven wide and two medium bands, notably including F480M, which enables clean isolation of the $H\alpha$ line at the highest redshift possible with NIRCam data.

The paper is structured as follows. Section 2 outlines our selection method, including dedicated simulations for completeness and contamination. In Section 3, we describe the observations and data processing of GLIMPSE and other public NIRCam data sets suitable for the Pop III galaxy search. In Section 4, we present Pop III candidates identified in the deep NIRCam data based on our new technique and discuss their physical properties and volume density. In Section 5, we explore an array of alternate scenarios regarding the origins of our Pop III galaxy candidates. We discuss future prospects in light of our findings in Section 6 and summarize our results in Section 7. Throughout this paper, we assume a flat Universe with $\Omega_m = 0.3$, $\Omega_\Lambda = 0.7$, $\sigma_8 = 0.8$, and $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$. We use magnitudes in the AB system (J. B. Oke & J. E. Gunn 1983). We adopt an uncertainty floor of 5% for the input photometry when running EAZY and SED fitting to capture systematic uncertainties (e.g., A. Weibel et al. 2024).

## 2. Efficient Pop III Search with NIRCam

In this section we establish that a Pop III search using NIRCam is feasible, and that such a search may be executed with high completeness and negligible contamination. We carefully quantify the imaging depth and filter-set coverage required for a successful search. We begin with Figure 1, where we compare the spectral energy distributions (SEDs) of Pop III and metal-enriched galaxies at $z=6.4$. Three key features of Pop III SEDs are highlighted: the absence of [O III], the strong $H\alpha$, and the Balmer jump.<fnref n="29" /> These three features are strong enough to appear in the NIRCam photometry out to $z \sim 7$ enabling us to develop an efficient selection method (see also A. K. Inoue 2011). We note that these SED features may diminish with age after the Pop III starburst or in case the escape fraction of ionizing photons is relatively high. Our study specifically focuses on enabling the search for the young phase of Pop III galaxies, rather than encompassing all of their evolutionary stages.

To quantitatively analyze these features and develop a search strategy, we introduce two photometric approaches purely using NIRCam in the following subsections: (i) color-color diagrams, and (ii) SED fitting. We also evaluate the completeness and contamination rates for each method, with the goal of identifying the most effective approach.

### 2.1. Color-based Selection

In Figure 2, we present possible color-color spaces for distinguishing Pop III galaxies from metal-enriched galaxies at $z \sim 6-7$. The filter combinations are optimized to focus on the

<footnote n="29">²⁹ Although the Balmer jump is not unique to Pop III stellar populations, it becomes more pronounced (e.g., [F200W - F410M] < -1.0; see Figure 2) in very young, metal-free systems due to their hotter stars. Therefore, in this study, we include the Balmer jump among the key features for distinguishing Pop III galaxies in broadband photometry.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Key Pop III galaxy features motivating our search strategy illustrated using the GLIMPSE filter-set. Top: Yggdrasil Pop III galaxy SED model (E. Zackrisson et al. 2011) with a moderately top-heavy IMF (blue curve). This model is characterized by $log(M_{star}/M_{\odot}) = 6$, $f_{cov} = 1.0$, and $t_{age} = 1.0$ Myr at $z = 6.4$, with a characteristic mass of individual Pop III stars of $M_{*,IMF} \simeq 10 M_{\odot}$ (Pop III.2). Middle: a metal-enriched galaxy SED model (green curve) is also shown, generated using BAGPIPES with the following parameters: $Z_{gas}/Z_{\odot} = 0.1$, $n_e = 1000$ cm$^{-3}$, and $log(U) = -2$. In both models, we masked the strong Ly$\alpha$ given its frequent damping by the neutral intergalactic medium (IGM) at this redshift. Bottom: the filter responses of the nine NIRCam filters used in GLIMPSE. In the top and middle panels, the expected NIRCam photometry is presented as open circles, where the y-axis errors denote the $1\sigma$ uncertainty based on the GLIMPSE data depth. Red shading highlights several unique SED features of the Pop III galaxy, including strong Balmer lines, the absence of the [O III] line, and the Balmer jump, that are effectively captured by the NIRCam photometry.</figcaption></figure>

distinctive characteristics of Pop III objects described above—the complete absence of [O III], strong hydrogen lines, and Balmer jump.

For Pop III galaxies, we utilize the Yggdrasil models presented in E. Zackrisson et al. (2011). The Yggdrasil framework combines stellar and nebular emission contributions by employing synthetic spectra from stellar population models and photoionization calculations. Additionally, we incorporate the Cloudy-based Pop III models introduced in K. Nakajima & R. Maiolino (2022), which explore a range of physical conditions for Pop III stellar populations and their surrounding gas. We summarize the key parameters of these models in Table 1 and provide their detail descriptions in Appendix A.

For modeling metal-enriched (albeit still quite metal-poor, “Population II” galaxies), we employ mock galaxy SEDs generated using BAGPIPES (A. C. Carnall et al. 2018), a Bayesian spectral fitting code designed for modeling galaxy emission across far-ultraviolet to microwave wavelengths, including stellar population synthesis models, nebular emission, and dust attenuation. The key parameters and detailed descriptions of the BAGPIPES-based mock SEDs are also summarized in Table 1 and Appendix A.

Figure 2 presents the color tracks of Pop III galaxies (blue symbols) at $z = 5.6–6.6$ and metal-enriched galaxies (other colored symbols) at $z = 5.0–7.5$. The wide redshift range displayed for the metal-enriched galaxies is helpful to ensure separation even with contaminants from different redshifts in the color–color spaces (see Section 5.3 for discussion of even

lower-redshift interlopers). We find that Pop III galaxies are distinctly separated from metal-enriched galaxies in specific regions of the color–color space owing to their unique spectral features highlighted in Figure 1.

It is important to note that these features, especially the strong H$\alpha$ line reflected in the y-axis colors in all diagrams, are sensitive to the stellar age ($t_{age}$). The colors of the Pop III galaxies suggest that these distinctive features begin to diminish around $t_{age} \simeq 10–30$ Myr (see faint cyan points in Figure 2). Importantly, the effectiveness of the color separation depends on the redshift range, as the technique relies on strong emission lines falling in specific NIRCam filters. Since Pop III star formation is more likely to occur at higher redshifts, this study focuses on the $z \simeq 6–7$ range, where the H$\alpha$ line is still captured by the NIRCam filters. However, extending this technique to lower redshifts using different filter sets is an avenue worth exploring in future studies (see $z \sim 4–5$; M. Nishigaki et al. 2023). In Appendix C, we also discuss the contamination rate when extending the Pop III search to lower redshifts. We also note that the separation between populations seen in Figure 2 is driven by F410M, which is sensitive to the continuum level, thereby constraining the high H$\alpha$ EW and Balmer jump, highlighting the crucial role played by medium bands. The separation could be made even cleaner when adding additional medium bands that help more securely constrain the line strengths of H$\alpha$ and [O III]. However, the effective redshift range probed becomes narrower to capture the specific emission lines with specific medium-band filters. In Appendix B, we discuss the color–color diagrams and the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** NIRCam color–color diagrams for selecting $z \simeq 6-7$ Pop III galaxies. In all diagrams, the selection (dashed red line) is guided by the same set of features: the x-axis colors detect weak [O III] lines and/or the Balmer jump, while the y-axis colors capture the strong H$\alpha$ EW. The cyan symbols represent Pop III galaxies at $z = 5.6-6.6$, taken from the Yggdrasil models (E. Zackrisson et al. 2011), assuming a very top-heavy initial mass function (IMF; Pop III.1), a moderately top-heavy IMF (Pop III.2), and a standard Kroupa IMF (Pop III.Kroupa) with $t_{age} = 0.01-30$ Myr. K. Nakajima & R. Maiolino (2022) models are also shown, assuming a Salpeter IMF with mass ranges of 1–100 M$_{\odot}$ (Pop III.Sal1-100), 50–500 M$_{\odot}$ (Pop III.Sal50-500), and 1–500 M$_{\odot}$ (Pop III.Sal1-500). The transparency of the symbols increases along three tiers of $t_{age}$: 0.01–10 Myr, 10–20 Myr, and 20–30 Myr. The other colored symbols represent metal-enriched mock galaxies at $z = 5.0-7.5$, generated using BAGPIPES. The blue, dark-magenta, pink, and gold colors denote $Z_{gas}/Z_{\odot} = 0.01, 0.05, 0.10$, and 0.20, respectively. The solid, dashed, and dotted lines correspond to $\log(U) = -1, -2$, and $-3$, respectively. All metal-enriched galaxies are assumed to have $A_{V} = 0$. The brown curve represents a dusty star-forming galaxy at $z = 5.0-7.5$ with $Z_{gas}/Z_{\odot} = 1.0, A_{V} = 2.0$, and $\log(U) = -2$ (its red color is off the scale in the middle panel). The gray curve indicates cool dwarfs with $T_{eff} = 500-2500$ K taken from Phoenix stellar template library in EAZY. In all symbols, the sizes, ranging from small to large, correspond to increasing redshift. Accounting for contaminants from a wide range of redshifts is critical in the selection process, and thus, this wide range is displayed for the metal-enriched galaxies.</figcaption></figure>

effective redshift range when utilizing F480M instead of F444W.

To facilitate the identification of Pop III candidates, we define the red dashed-line regions in the color–color diagrams as selection criteria for Pop III galaxies. Since this color selection is effective only within specific redshift ranges, we first perform an SED template fitting using EAZY (G. B. Brammer et al. 2008) with the 14 Pop III templates described in Section 2.2. We note that this means these two steps of redshift and color–color selections are not strictly independent, but this is essential to ensure the applicability of the color diagrams. We set the following redshift-based selection criteria:

$$5.0 < z_{\text{phot,PopIII}} < 7.5 \land z_{\text{phot,PopIII}} - z_{16} < 0.5 \land z_{84} - z_{\text{phot,PopIII}} < 0.5, \quad (1)$$

where $z_{\text{phot,PopIII}}$ refers to the best-fit photometric redshift derived from EAZY with the Pop III templates, while $z_{16}$ and $z_{84}$ denote the redshift at the 16th and 84th percentiles. Then, we step forward to the color selection enclosed by the following vertices in the color–color space:

1. F277W – F356W versus F410M – F444W (vertices):

$$(0.5, 2.5), (0.3, 1.0), (-0.12, 0.9), (-0.12, 0.0), (-0.8, 0.0), (-0.8, 2.5) \quad (2)$$

2. F200W – F410M versus F356W – F444W (vertices):

$$(0.6, 1.2), (0.4, 0.65), (-0.55, 0.55), (-0.8, -0.1), (-2.5, -0.1), (-2.5, 1.2) \quad (3)$$

3. F356W – F410M versus F410M – F444W (vertices):

$$(1.5, 0.0), (0.6, 0.0), (-0.2, 0.8), (-2.5, 2.0), (-2.5, 2.7), (1.5, 2.7), (1.5, 1.8). \quad (4)$$

These three color diagrams are based on the same principles: the x-axis colors are sensitive to the weak [O III] line and/or the Balmer jump, while the y-axis colors are designed to detect the strong H$\alpha$ line (see also Figure 1). Due to the clear separation observed in all three color–color diagrams, any one of these criteria alone could suffice for identifying Pop III candidates. However, applying all three criteria together helps mitigate contamination from noise fluctuations, especially for low signal-to-noise ratio (S/N) sources. Each diagram places different weights on specific features: for instance, the [F200W – F410M] color in the second diagram is more directly sensitive to the Balmer jump, the [F356W – F410M] color in the third diagram better captures the weak [O III] line, and the [F277W – F356W] color in the first diagram is sensitive to both features. Therefore, combining these three criteria enhances the robustness of the selection process.<fnref n="30" /> Next, we discuss the SED-based selection method after which we turn to evaluate the completeness (Section 2.3) and contamination rate (Section 2.4) to ensure reliable selection of Pop III candidates.

<footnote n="30">³⁰ In the future, the application of machine learning techniques to perform multidimensional separation could further improve the selection process (e.g., T. Kojima et al. 2020).</footnote>
<table><caption>**Table 1**
Summary of the Models Used in Our Selection Methods, Completeness/Contamination Simulations, and Additional Tests</caption><thead><tr><th>Type</th><th>Reference</th><th>IMF</th><th>Key Descriptions</th><th>Usage</th></tr></thead><tbody><tr><td>Pop III</td><td>E. Zackrisson et al. (2011)<br>K. Nakajima & R. Maiolino (2022)</td><td>Power-law (α = 2.35) 50–500M<sub>☉</sub> (PopIII.1), Log-normal (center = 10M<sub>☉</sub>, σ = 1M<sub>☉</sub>) 1–500M<sub>☉</sub> (PopIII.2), Kroupa 0.1–100M<sub>☉</sub> (PopIII.Kroupa)<br>Salpeter 1–100M<sub>☉</sub> (PopIII.Sal1-100), 1–500M<sub>☉</sub> (PopIII.Sal1-500), 50–500M<sub>☉</sub> (PopIII.Sal50-500)</td><td>n<sub>e</sub> = 100 cm<sup>-3</sup>, f<sub>cov</sub> = 1.0, t<sub>age</sub> = [0.01, 0.1, 1.0, 3.6] Myr (PopIII.1), [0.01, 1.0, 5, 10, 50] Myr (PopIII.2), [0.01, 1.0, 5, 10, 50] Myr (PopIII.Kroupa)<br>n<sub>e</sub> = 10<sup>3</sup> cm<sup>-3</sup>, log U = -1.5, t<sub>age</sub> = 0 Myr (zero age)</td><td>Color-based selection (Section 2.1), SED-based selection (Section 2.2), Completeness simulation (Section 2.3)<br>SED-based selection (Section 2.2)</td></tr><tr><td>Pop II</td><td>BAGPIPES (A. C. Carnall et al. 2018)<br>EAZY (G. B. Brammer et al. 2008)<br>SC SAM (L. Y. A. Yung et al. 2022)</td><td>Kroupa<br>Salpeter<br>Chabrier</td><td>Z<sub>gas</sub>/Z<sub>☉</sub> = [0.01, 0.05, 0.10, 0.20] log U = [-3, -2, -1],<br>A<sub>V</sub> = 0<br>blue_sfhz_13, a strong line emitter template (A. C. Carnall et al. 2023b)<br>Semianalytically simulated 1,472,791 galaxies across a light cone in ~800 arcmin<sup>2</sup> at z = 0–10</td><td>Color-based selection (Section 2.1)<br>SED-based selection (Section 2.2)<br>Contamination simulation (Section 2.2)</td></tr><tr><td>Metal-poor AGN</td><td>K. Inayoshi et al. (2022a, 2022b)<br>K. Nakajima & R. Maiolino (2022)</td><td>...</td><td>Seed BH mass M<sub>BH</sub> = 10<sup>6</sup> M<sub>☉</sub> in a metal-poor galaxy (Z<sub>gas</sub>/Z<sub>☉</sub> = 0.01). Broken power-law accretion disk, nebular emission, and dense accretion disk radiation (0.1–100 pc), with viewing angle 60°<br>DCBH with thermal Big Bump (T<sub>bb</sub> = 5 × 10<sup>4</sup>, 1 × 10<sup>5</sup>, 2 × 10<sup>5</sup> K), power-law continuum (α = -1.2, -1.6, -2.0), and nebular emission (n<sub>e</sub> = 10<sup>3</sup> cm<sup>-3</sup>, Z<sub>gas</sub>/Z<sub>☉</sub> = 0, log U = -1.5) in plane parallel</td><td>Additional SED test (Section 5.2, Appendix F)<br>Additional SED test (Section 5.2, Appendix F)</td></tr></tbody></table>
## 2.2. SED-based Selection

Similar to the high-redshift galaxy search (e.g., S. L. Finkelstein et al. 2015), another promising method for identifying Pop III candidates is using SED fitting, in addition to the color-based approach (Section 2.1). This SED-based approach leverages the full photometric data set rather than relying on specific filter combinations optimized for certain redshift ranges. By maximizing the information extracted from all available photometric bands, SED fitting may enable a more comprehensive search for robust Pop III candidates. We extend the concept of robustness of high-redshift candidates based on the difference in the $\chi^2$ between the best-fit high-z and forced low-z SEDs to the Pop III use case (e.g., R. P. Naidu et al. 2022b; C. T. Donnan et al. 2023; Y. Harikane et al. 2023a; S. L. Finkelstein et al. 2024). We use EAZY (G. B. Brammer et al. 2008) to compute the $\chi^2$ difference between the best-fit metal-enriched galaxy templates and Pop III templates.

The first template set consists of the blue SFH 13 model subset,<fnref n="31" /> aiming to obtain the best-fit SED as a metal-enriched galaxy. This set incorporates redshift-dependent star formation histories (SFHs) and dust attenuation. Additionally, the linear combinations of log-normal SFHs are further constrained to ensure they do not start earlier than the age of the Universe. The template set is further expanded with a blue galaxy template derived from a high-SNR JWST spectrum of a lensed $z=8.50$ galaxy representative of typical nebular conditions in the early Universe (A. C. Carnall et al. 2023a). This template extends the color grid coverage for galaxies with extreme emission-line equivalent widths, resulting in a total of 14 templates.

The second template set comprises 14 Pop III templates based on the Yggdrasil model (E. Zackrisson et al. 2011), aiming to identify the best-fit SED of a Pop III galaxy. These templates are selected from three IMF scenarios of Pop III.1, Pop III.2, and Pop III.Kroupa summarized in Table 1. From Pop III.1 (the most top-heavy IMF), we use three templates with stellar ages $t_{age} = [0.01, 1.0, 3.6]$ Myr. Pop III.2 (moderately top-heavy IMF) contributes five templates with $t_{age} = [0.01, 1.0, 5, 10, 50]$ Myr, and Pop III.Kroupa (Kroupa IMF) provides six templates with similar stellar ages of $t_{age} = [0.01, 1.0, 5, 10, 50, 100]$ Myr. Note that the difference in age range across the IMF flavors is because no massive stars survive beyond a few Myr in the top-heavy scenario. In all cases, the gas covering fraction is fixed at $f_{cov} = 1.0$, ensuring maximal nebular contribution. Although additional Pop III templates, such as those from K. Nakajima & R. Maiolino (2022), could be included, their contribution to the color-space coverage is minimal (see Figure 2). To maintain balance with the first template set, we limit the second set to 14 templates, ensuring the same number of models for a fair evaluation of the delta $\chi^2$ measurement. While we acknowledge that additional Pop III templates are available in the Yggdrasil library, we find that increasing the number of templates in EAZY does not necessarily improve photometric redshift performance, and can in fact introduce greater scatter. We thus adopt this moderate set of 14 templates, which spans a representative range of Pop III models with different IMFs and ages (0.01–100 Myr), and consider finer sampling of the

Pop III parameter space to be beyond the scope of this initial systematic NIRCam search.

With these template sets, we run EAZY and obtain the $\chi^2$ values for the best-fit metal-enriched galaxy template, $\chi^2_{galaxy}$, and Pop III galaxy template set, $\chi^2_{PopIII}$. By defining $\Delta\chi^2 \equiv \chi^2_{galaxy} - \chi^2_{PopIII}$, we set the following $\chi^2$ criteria:

$$ \Delta\chi^2 \ge C_{thresh} \wedge \chi^2(\text{PopIII}) < 10 \quad (5) $$

or

$$ \Delta\chi^2 \ge 30 \wedge \chi^2(\text{PopIII}) < 20, \quad (6) $$

where $C_{thresh}$ is a threshold value. In high-redshift galaxy studies (e.g., R. P. Naidu et al. 2022b; C. T. Donnan et al. 2023; Y. Harikane et al. 2023a; S. L. Finkelstein et al. 2024), $C_{thresh}$ is typically set to either 4 or 9.<fnref n="32" /> The choice of $C_{thresh}$ represents a trade-off between completeness and contamination rates. Given the scope of our search for extremely rare objects like Pop III galaxies, we prioritize mitigating the contamination rate and adopt $C_{thresh} = 9$ in the following analyses. Since the $\chi^2$ value can naturally increase with the number of filters used in SED fitting, Equation (6) is incorporated to ensure that potential candidates are not overlooked in regions where more extensive filter sets are available. Unlike the color-based selection (Section 2.1), the SED-based selection is not restricted to a specific redshift range. However, to ensure a fair comparison between the two methods, we adopt the same redshift criteria defined in Equation (1). In Appendix C, we discuss the distributions of the completeness and contamination rate when the SED-based method is applied to a broader redshift range.

Note that we use the photometric redshift estimated with the Pop III templates ($z_{phot,PopIII}$), rather than galaxy templates, for the redshift criteria above. This choice is motivated by the fact that the Lyman-break feature is not always well constrained at this redshift range, particularly in the faint regime (>28 mag), where extremely deep HST/Advanced Camera for Surveys (ACS) data would be required. Without robust constraints from the Lyman-break feature, general galaxy templates may produce high $\chi^2$ values due to their inability to reproduce the unique Pop III signatures (Figure 1), potentially leading to incorrect redshift estimates. While Pop III templates may similarly result in incorrect redshift estimates for general metal-enriched galaxies, the primary goal of this analysis is to identify Pop III candidates, rather than providing the most accurate redshift estimates for the majority of metal-enriched galaxies. Therefore, the photometric redshift refers to the estimate based on the Pop III templates throughout this paper (otherwise specified), including the following completeness and contamination rate measurements.

## 2.3. Completeness

To evaluate the selection function of our color-based and SED-based methods, we assess the completeness of each approach using Monte Carlo (MC) simulations. For both simulations, we adopt the Yggdrasil Pop III.2 model

<footnote n="31">³¹ https://github.com/gbrammer/EAZY-photoz/tree/master/templates/sfhz</footnote>

<footnote n="32">³² These thresholds correspond approximately to $2\sigma$ and $3\sigma$ significance levels for a degree of freedom (dof) of 1. However, as SED fitting involves multiple parameters beyond redshift, the actual dof is >1, and this statistical interpretation should be regarded as an approximation.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Completeness (left) and contamination rate (right) of the Pop III selection based on the color-based (top; Section 2.1) and SED-based (bottom; Section 2.2) methods. The cyan contours denote 0.5 and 0.8 for completeness and 0.01, 0.1, and 0.5 for contamination rate. For completeness, we use a Yggdrasil Pop III model (Pop III.2, $t_{age}$ = 1.0 Myr, $f_{cov}$ = 1.0), while confirming similar results with other templates representing young stellar populations (i.e., those with strong nebular emission lines). For contamination rate, we utilize the Santa Cruz Semi-Analytic Model (SC SAM) catalog (L. Y. A. Yung et al. 2022), which includes 1,472,791 metal-enriched mock galaxies spanning $z$ = 0–10 in a light cone covering ≃800 arcmin². To ensure a fair assessment, we first derive $z_{phot,PopIII}$ using EAZY and bin the galaxies in redshift and observed F200W magnitude grids. The completeness and contamination rate are then calculated in each grid through Monte Carlo simulations by adding random noise based on the GLIMPSE data depth and applying the same selection procedures used in the observations. The completeness and contamination rates reflect a trade-off between the two methods. The color-based method achieves high completeness (>50%–80%) within a limited redshift range, accompanied by a low contamination rate, whereas the SED-based method provides high completeness over a broader redshift range but with a relatively higher contamination rate. While both completeness and contamination rates are critical, we prioritize achieving a low contamination rate in the search for extremely rare objects like Pop III galaxies. To this end, we combine the color-based and SED-based methods to maximally mitigate contamination while maintaining robust completeness (Figure 4).</figcaption></figure>

(E. Zackrisson et al. 2011) with $t_{age}$ = 1.0 Myr and $f_{cov}$ = 1.0. Note that we verify negligible differences in the final completeness results when repeating the simulations with various Pop III models, provided $t_{age} \lesssim$ 15 Myr. This outcome is expected, as the key SED features (see Figure 1) begin to diminish as $t_{age}$ increases (see Figure 2). NIRCam photometry is calculated for the nine filters used in GLIMPSE, shifting and rescaling the template across redshift and $m_{F200W}$ grids spanning $z$ = [5.0:7.5] and $m_{F200W}$ = [31.5:28.0], respectively. At each grid point, 100 mock NIRCam photometric sets are generated by injecting random noise based on the detection limits of each filter. To ensure consistency in redshift criteria across both methods, we apply the $z_{phot}$ selection defined in Equation (1) by running EAZY with the Pop III template set for the 100 mock photometric sets at each grid point. For the SED-based method, we also separately run EAZY with the metal-

enriched galaxy templates. Following the selection procedures described in Sections 2.1 and 2.2, we calculate the completeness for each grid point for both methods.

The left panels of Figure 3 show the completeness estimates for the color-based (top panels) and SED-based (bottom panels) methods. For the color-based method, we find that the completeness function exhibits steep cutoffs at $z \simeq 5.6$ and $z \simeq 6.6$. This reflects the dependence of the color-based method on the specific NIRCam filters to capture the unique SED features of Pop III galaxies, such as bright hydrogen lines and the absence of the [O III] line (Figure 1). In other words, at $z \lesssim 5.6$ or $z \gtrsim 6.6$, either the [O III] or H$\alpha$ emission line falls outside the F356W and F444W filters, resulting in a loss of these diagnostic features. Within the $z \simeq 5.6–6.6$ range, the completeness exceeds ~50%–80% for $m_{F200W} \le 30.0$, ensuring high completeness within this specific redshift range. In
contrast, the SED-based method shows high completeness over a broader redshift range, reflecting its ability to leverage the full photometric data set rather than being restricted to specific filters optimized for a particular redshift. However, the completeness decreases at the faint end ($m_{F200W} > 30.0$ mag) for $z \gtrsim 6.6$, as the Hα line moves out of the F444W filter, losing one of the strongest features of Pop III galaxies in the SED.

## 2.4. Contamination Rate

Here we investigate the contamination rate for both methods. For this analysis, we utilize the DR3 mock galaxy catalog generated with the Santa Cruz Semi-Analytic Model, as presented in L. Y. A. Yung et al. (2022).<fnref n="33" /> This catalog contains 1,472,791 metal-enriched galaxies across a light cone spanning an area of ~800 arcmin² and covering the redshift range $z=0–10$. The catalog includes photometry for all NIRCam filters, accounting for emission-line contributions from the mock galaxy spectra. To simulate realistic observations, we inject random noise into the NIRCam photometry based on the GLIMPSE detection limits (see Section 3.1). We then derive the observed $m_{F200W}$ magnitudes and estimate $z_{phot,PopIII}$ using EAZY. The mock galaxies are sorted into redshift and $m_{F200W}$ grids spanning $z = [5.0:7.5]$ and $m_{F200W} = [31.5:28.0]$. Following the selection criteria defined in Sections 2.1 and 2.2, we estimate the contamination rate for each grid in both methods. Note that we confirm that most grids contain >10–100 galaxies, ensuring reliable statistics.

The right panels of Figure 3 show the contamination rate—the fraction of metal-enriched galaxies that satisfy the Pop III criteria—estimates for the color-based (top panels) and SED-based (bottom panels) methods. For the color-based method, the contamination rate is nearly zero across most of the parameter space, with only a few regions at the faintest regime ($\simeq 30–31$ mag) exhibiting modest contamination rates ($\sim 1\%$). This is likely due to the noise fluctuation at the faintest regime close to the detection limit, while this modest contamination rate highlights the robustness of the color-based method. In contrast, the SED-based method displays contamination rates exceeding $\sim 1\%$ across almost the entire parameter space, with rates exceeding 10%–50% in certain regions.<fnref n="34" /> However, within the specific redshift range $z \simeq 5.8–6.6$, the contamination rate remains modest ($< 10\%$). Similar to the color-based method, these are the redshifts where the unique SED features of Pop III galaxies—strong Balmer lines and the absence of the [O III] line—are effectively captured by the NIRCam filters.

Although the SED-based method achieves high completeness over a broad redshift range, there is significant contamination from metal-enriched galaxies. To mitigate this, it is necessary to focus on an optimal redshift range that is very similar to the color-based method. This motivates our choice to apply the SED-based method to decrease contamination by adding it as an additional selection criterion to the color-based method, rather than applying it independently across a wide redshift range.

## 2.5. Fiducial Selection Method

Based on the results of the completeness and contamination tests, we combine the selection criteria of the color-based and SED-based methods (i.e., Pop III candidates must satisfy Equations (1)–(6)). Figure 4 shows the completeness and contamination rates derived from the combined color+SED method. The completeness remains largely similar to that of the color-based method due to the broader selection range allowed by the SED-based method. However, the contamination rate is reduced to nearly zero across almost the entire parameter space, achieved through the strict selection function provided by the combination of both methods. We therefore adopt this combined technique as our fiducial selection method.

## 2.6. Required Depth for the Pop III Search

Although we define the fiducial selection function in Section 2.5, the selection focuses on the SED colors alone. It is also important to consider the possible luminosity for the Pop III galaxies and understand the required data depth for a meaningful Pop III search. The first galaxies are predicted to form within halos where the cooling of primordial gas enables star formation. At $z \sim 6–7$, the increasing background of Lyman–Werner radiation likely suppresses molecular hydrogen cooling, making atomic cooling the dominant mechanism. In these cases, a virial temperature of $T_{vir} \sim 10^4$ K, the threshold for efficient atomic cooling, corresponds to a halo mass of $M_{vir} \sim 10^8 M_{\odot}$ (e.g., V. Bromm & N. Yoshida 2011; M. Kulkarni et al. 2021). In regions yet to be reionized, such metal-free halos are expected to form Pop III stars upon reaching the atomic cooling limit. Assuming a nominal gas conversion factor $f_{cool} = 0.01$ within halos and a time-averaged star formation efficiency $\epsilon_{*} = 0.1$, a single starburst is expected to yield a stellar mass of $M_{star} \simeq 10^5 M_{\odot}$ (see, e.g., review by V. Bromm & N. Yoshida 2011). This estimate may be conservative, as recent JWST studies at $z \gtrsim 6$ suggest $\epsilon_{*}$ values $> 0.1$. This is supported by the significant abundance of UV-bright galaxies at $z \gtrsim 9$ (e.g., Y. Harikane et al. 2023a; S. L. Finkelstein et al. 2024), the early formation of massive galaxies by $z = 5–6$ (e.g., A. de Graaff et al. 2025; A. Weibel et al. 2024; M. Xiao et al. 2024), and the high stellar and gas densities observed in resolved lensed galaxies at $z = 6–10$ (e.g., A. Adamo et al. 2024; S. Fujimoto et al. 2024a; M. Dessauges-Zavadsky et al. 2025; however, see also C. T. Donnan et al. 2025, for an alternative viewpoint). Nevertheless, our simple estimate provides a useful baseline.

Figure 5 illustrates how the apparent magnitude of Pop III galaxies varies as a function of $t_{age}$. We use the Yggdrasil Pop III galaxy template ($M_{star} = 10^6 M_{\odot}$) with three different IMFs, and also present cases for $10^5 M_{\odot}$ and $10^4 M_{\odot}$ by assuming simple scaling, as per the baseline $M_{star}$ for first galaxies that are dominated by Pop III stars as described above. For reference, the typical $5\sigma$ depths of JWST legacy surveys are shown with horizontal gray lines, corresponding to three survey layers: wide ($\simeq 28$ mag), medium-deep ($\simeq 29$ mag), and deep ($\simeq 30.5$ mag) surveys (e.g., C. M. Casey et al. 2023; M. B. Bagley et al. 2023; C. T. Donnan et al. 2024; S. L. Finkelstein et al. 2024; B. Robertson et al. 2024). In most cases, the Pop III galaxies are likely to be too faint ($\gtrsim 30.5$ mag) for these JWST surveys, except under the most optimistic conditions, such as $M_{star} \simeq 10^6 M_{\odot}$, top-heavy

<footnote n="33">³³ The catalog is available at https://ceers.github.io/sdr3.html#catalogs.</footnote>

<footnote n="34">³⁴ In the contamination rate distribution shown in Figure 3, we note that the redshift is binned based on $z_{phot,PopIII}$, rather than the true redshift of the mock galaxies, to ensure consistency with the procedure applied to the observed data. Therefore, the high contamination rates of > 10%–50% in specific parameter spaces only apply to sources whose SEDs are initially favored by the Pop III templates at $z \simeq 5.0–7.5$.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Completeness and contamination rate for our fiducial Pop III selection combining the color-based and SED-based methods. The contamination rate is reduced to zero across the entire parameter space, ensuring a robust search for Pop III galaxies with NIRCam. The red pentagon represents the parameter space of the promising Pop III candidate GLIMPSE-16043 found in GLIMPSE that satisfies our fiducial selection criteria.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** A meaningful Pop III search demands the deepest JWST surveys as well as the assistance of gravitational lensing. Apparent magnitudes in the rest-frame UV (solid lines; F150W) and optical (dashed lines; F410M) are shown as a function of $t_{age}$ for Pop III galaxies at $z = 6.5$. The three colored curves represent the Yggdrasil Pop III models with different IMFs: Pop III.1 (blue; very top-heavy), Pop III.2 (orange; moderately top-heavy), and Pop III.Kroupa (green; standard). The three $M_{star}$ cases are shown from left to right, with a simple scaling applied to the Yggdrasil Pop III models, which originally assume $M_{star} = 10^6 M_{\odot}$. Following the arguments outlined in Section 2.6, we consider $\approx 10^5 M_{\odot}$ as our fiducial case—this makes clear that only the deepest JWST surveys, further aided by lensing, have a reasonable chance at detecting Pop III stars. While the mass-to-light ratio could vary with $M_{star}$ due to differences in nebular parameters (e.g., electron density, ionization parameter), this figure provides a first-order insight into the NIRCam depth required for a meaningful Pop III search (see Section 2.6). The future prospects for the meaningful survey design are also discussed in Section 6.</figcaption></figure>

IMFs, and $t_{age} \lesssim 5$ Myr. This highlights the immense challenge in successfully identifying Pop III galaxies, consistent with the absence of robustly confirmed Pop III sources to date. From Figure 5, it is clear that leveraging gravitational lensing is critical to push past $\approx 30.5$ mag, without which the majority of scenarios remain out of reach (involving, e.g., slightly older stellar ages $> 5$ Myr). Indeed, the most metal-poor source detected with JWST until date is a highly magnified clump at $z = 6.6$ ($\mu > 100$; E. Vanzella et al. 2023a).

### 3. Application to Deep JWST Surveys

In Section 2, we outlined our novel NIRCam selection method for selecting $z \simeq 6-7$ Pop III galaxies and the required depths for a meaningful search. Based on the filters used for

the color-color diagrams, which are optimized to capture the key features of the Balmer jump, strong H$\alpha$ emission, and the absence of the [O III] line, the NIRCam selection method requires the F200W, F277W, F356W, F410M, and F444W filters (Figure 2). Also given the preference for deeper data (Figure 5), we focus on the following publicly available NIRCam data from JWST legacy surveys: CEERS (S. L. Finkelstein et al. 2024), PRIMER-UDS/COSMOS (C. T. Donnan et al. 2024), JADES Origins Field (JOF; B. Robertson et al. 2024), UNCOVER+MegaScience (R. Bezanson et al. 2024; K. A. Suess et al. 2024), and the latest Cycle 2 large program, GLIMPSE (H. Atek et al. 2025, in preparation). We summarize the data sets used in our analysis in Table 2. Collectively, these surveys are representative of the various depth tiers shown in Figure 5, enabling a systematic “wedding
<table><caption>**Table 2**
Summary of Data Sets Used in Our Systematic Pop III Galaxy Search</caption><thead><tr><th>Field</th><th>NIRCam Surveys</th><th>Area<br>(arcmin²)</th><th>Depth<sup>a</sup><br>(AB mag)</th><th>Filterset<sup>b</sup></th></tr></thead><tbody><tr><td>Abell S1063</td><td>GLIMPSE (H. Atek et al. 2025, in preparation; this work; V. Kokorev et al. 2025)</td><td>9.4</td><td>30.2–30.5<sup>c</sup></td><td>≥0.9 µm wide filters + F410M + F480M</td></tr><tr><td>A2744</td><td>UNCOVER (R. Bezanson et al. 2024), MegaScience (K. A. Suess et al. 2024), ALT (R. P. Naidu et al. 2024), BEACONS (T. Morishita et al. 2024), #2883 (PI: Sun), #2756 (PI: Chen), #3538 (PI: Iani)</td><td>45</td><td>29.2–30.2<sup>c</sup></td><td>All NIRCam wide and medium bands</td></tr><tr><td>GOODS-S</td><td>JOF (B. Robertson et al. 2024), JADES (D. J. Eisenstein et al. 2023a, 2023b), FRESCO (P. A. Oesch et al. 2023)</td><td>9.4</td><td>30.2–31.0</td><td>≥0.9 µm wide filters + F182M + F210M + F250M + F300M + F335M + F410M</td></tr><tr><td>EGS</td><td>CEERS (S. L. Finkelstein et al. 2024)</td><td>88.1</td><td>28.6–29.2</td><td>≥1.1 µm wide filters + F410M</td></tr><tr><td>UDS</td><td>PRIMER (C. T. Donnan et al. 2024)</td><td>170.7</td><td>27.6–28.4</td><td>≥0.9 µm wide filters + F410M</td></tr><tr><td>COSMOS</td><td>PRIMER (C. T. Donnan et al. 2024)</td><td>132.2</td><td>27.6–28.4</td><td>≥0.9 µm wide filters + F410M</td></tr></tbody></table>

**Notes.**

cake” search. Below, we briefly describe the various surveys and how we incorporate their data.

### 3.1. GLIMPSE

We performed our own reductions for the GLIMPSE data, and here we briefly explain the survey, observations, data reduction, and its processing. The GLIMPSE survey (PID 3293; PIs: H. Atek & J. Chisholm) is designed to probe the faintest galaxies during the epoch of reionization. Observations were conducted between 2024 September 20 and September 28. The target was the well-studied massive lensing cluster Abell S1063, especially taking advantage of the largest number of highly magnified sources ($\mu > 10$), leveraged by one of the most robust lens models, among the six Hubble Frontier Field Clusters (e.g., J. M. Lotz et al. 2017). The program obtained ultradeep NIRCam imaging across seven broad bands at $\ge 0.9$ µm and two medium bands of F410M and F480M, with exposure times for each filter ranging from 16–39 hr. Specifically, the deepest imaging was performed in the F090W and F115W bands, with total exposure times of 39 hr each, while other filters had exposure times of ~16–23 hr, achieving the almost homogeneous sensitivity at ~1–5 µm. The availability of deep HST/ACS imaging (5σ depth of ~29 mag) is also beneficial for our study, providing crucial constraints on the Lyman-break features for our targeted populations at $z \sim 6–7$. The orientation of the observations was set to ensure optimal coverage of the lensing cluster. Module B was centered on the high-magnification core of the lensing cluster (magnification factor $\mu \gtrsim 3$), while Module B still covered moderately magnified regions ($\mu \sim 1.2–1.5$). GLIMPSE reaches depths of 30.2–30.5 AB mag (5σ, point source, 0.′2-diameter aperture) in broad bands, comparable to the deepest blank-field surveys such as NGDEEP (e.g., M. B. Bagley et al. 2023) and JOF (e.g., B. Robertson et al. 2024). By leveraging the magnification provided by Abell S1063, the effective sensitivities achieve unprecedented depths of $\gtrsim 33$ AB mag in the most highly

magnified regions ($\mu \gtrsim 10$), enabling the detection of galaxies an order-of-magnitude fainter than previously achievable.

The raw NIRCam imaging data were processed using the JWST pipeline (v1.9.7) with the 1293.pmap context map. Key steps included custom corrections for cosmic rays, 1/f noise, and detector artifacts, as well as refined background subtraction to mitigate contamination from bright cluster galaxies. The final mosaics were created using a drizzling algorithm, yielding pixel scales of 0.′02 for the short-wavelength channels and 0.′04 for the long-wavelength channels.

To account for contamination from intracluster light (ICL) and bright cluster galaxies (BCGs), we employed an iterative modeling approach. The ICL was modeled using smooth polynomial functions, while the BCGs were subtracted using GALFIT, employing Sérsic profile fits. These corrections were applied consistently across all exposures and filters. After removing ICL and BCG contributions, a secondary local background subtraction was performed to refine the photometry in affected regions. This multistep approach ensures accurate measurements of faint sources, particularly near the lensing cluster.

Sources were identified using SExtractor, with detections based on a combined image stack of the F277W, F356W, and F444W bands. To ensure consistent photometry across all bands, point-spread functions (PSFs) were empirically constructed by stacking isolated stars in the field. The PSFs of all bands were then matched to the lowest-resolution filter, F480M (FWHM $\simeq 0.′16$), using convolution kernels. Fluxes were measured in circular apertures with diameters ranging from 0.′1 to 1.′2, and aperture corrections were applied based on the empirically derived growth curves of the matched PSFs. We adopt the 0.′2-diameter aperture-corrected photometry in the following analyses. Photometric uncertainties were estimated by placing random apertures in source-free regions of the mosaics, thereby accounting for correlated noise and background variations. The final photometric catalog includes all sources with an S/N > 3 in at least three bands. Further

<footnote><sup>a</sup> All depths quoted are 5σ depths measured in 0.′2 diameter, with the exception of A2744 (0.′16 and 0.′32 diameter apertures for the SW and LW; R. Bezanson et al. 2024; K. A. Suess et al. 2024).</footnote>

<footnote><sup>b</sup> By selection, coverage in almost all SW+LW wide filters (F090W, F115W, F150W, F200W, F277W, F356W, and F444W) and at least one medium band (F410M) is available in all of our studied fields. A2744 is additionally covered by deep F070W imaging as well as by every single medium band (F140M, F162M, F182M, F210M, F250M, F300M, F335M, F360M, F410M, F430M, F460M, and F480M).</footnote>

<footnote><sup>c</sup> In fields leveraging gravitational lensing (Abell S1063, A2744), intrinsically fainter sources are accessible—e.g., the effective depth of GLIMPSE in Abell S1063 reaches >33 mag in regions with $\mu > 10$.</footnote>
details of the observation setup, scope, and data processing will be presented in a separate paper (H. Atek et al. 2025, in preparation; see also V. Kokorev et al. 2024).

## 3.2. UNCOVER+MegaScience

The UNCOVER survey (PID 2561; PIs: I. Labbé & R. Bezanson; R. Bezanson et al. 2024) and the “Medium Bands, Mega Science” follow-up program (PID 4111; PI: K. Suess; K. A. Suess et al. 2024) provide deep NIRCam imaging over the A2744 galaxy cluster, complementing earlier HST observations from the Frontier Fields program (J. M. Lotz et al. 2017). We utilize the publicly available UNCOVER DR3 (K. A. Suess et al. 2024) that additionally incorporates imaging from various NIRCam surveys including ALT (R. P. Naidu et al. 2024), BEACONS (T. Morishita et al. 2024), MAGNIF (#2833, PI: Sun), #2756 (PI: Chen), and #3538 (PI: Iani). The DR3 release delivers imaging in all 12 medium bands and all eight broadband filters, covering a continuous wavelength range from 0.7–5 µm. The depths span 29.2–30.2 mag (5σ, point source) for broad bands with 0.′′08-radius and 0.′′16-radius apertures for SW and LW filters, respectively. The inclusion of medium-band filters allows for precise sampling of key spectral features, such as the Balmer break, strong nebular emission lines (e.g., [O III] and Hα), and the rest-frame optical continuum at our redshifts of interest. This data set crucially benefits from lensing magnification (mean µ ≈ 2.5; L. J. Furtak et al. 2023; S. H. Price et al. 2024) provided by the A2744 cluster. The photometric catalog, reduced mosaics, and associated materials are publicly available from the UNCOVER and MegaScience project website.<fnref n="35" />

## 3.3. CEERS, PRIMER, and JOF

For the CEERS, PRIMER, and JOF regions, we utilize the photometric catalogs constructed in A. Weibel et al. (2024). The reduced images for these data sets are based on the DAWN JWST Archive’s v7 imaging products,<fnref n="36" /> which were processed using a custom JWST pipeline with standard calibrations, cosmic-ray removal, and background subtraction, ensuring consistent quality across fields (see more details in F. Valentino et al. 2023). The CEERS program (PID 1345; S. L. Finkelstein et al. 2024) offers imaging across the Extended Groth Strip (EGS) field, covering ~100 arcmin². Observations were conducted in seven NIRCam bands at ≥1.1 µm, including F410M, reaching depths of 28.6–29.2 mag (5σ, point source, 0.′′2-diameter; S. L. Finkelstein et al. 2024). The PRIMER survey (PID 1837; C. T. Donnan et al. 2024) is conducted in the COSMOS and UDS fields, with imaging over ~300 arcmin². Eight NIRCam filters at ≥0.9 µm are used, including F410M, with median depths of 27.6–28.4 mag (5σ, point source, 0.′′3-diameter; C. T. Donnan et al. 2024). JOF (PID 1180; B. Robertson et al. 2024) targets the GOODS-S region, offering ultradeep NIRCam imaging in 13 filters at ≥0.9 µm. The JOF data set achieves depths of up to 30.2–31.0 AB magnitudes (5σ, point source, 0.′′2-diameter; B. Robertson et al. 2024).

The photometry across these fields was systematically extracted. Source detection was performed using a weighted

<table><caption>**Table 3**
NIRCam Photometry of GLIMPSE-16043</caption><tr><td>R.A. (deg)</td><td>342.2123718</td></tr><tr><td>Decl. (deg)</td><td>-44.528751</td></tr><tr><td>F090W (nJy)</td><td>3.17 ± 0.60</td></tr><tr><td>F115W (nJy)</td><td>4.44 ± 0.52</td></tr><tr><td>F150W (nJy)</td><td>4.74 ± 0.61</td></tr><tr><td>F200W (nJy)</td><td>3.54 ± 0.62</td></tr><tr><td>F277W (nJy)</td><td>3.29 ± 0.71</td></tr><tr><td>F356W (nJy)</td><td>3.94 ± 0.74</td></tr><tr><td>F410M (nJy)</td><td>0.71 ± 1.21</td></tr><tr><td>F444W (nJy)</td><td>4.85 ± 0.67</td></tr><tr><td>F480M (nJy)</td><td>12.54 ± 2.47</td></tr></table>

**Note.** The measurements above are aperture-corrected photometry derived from PSF-matched images in the observed frame (i.e., without correcting for magnification).

combination of F277W, F356W, and F444W images. PSFs were constructed empirically and matched across all bands, enabling precise photometric measurements. Given the presence of massive, large galaxies in wide-area general fields, total fluxes were measured within 0.′′32-diameter apertures, which is slightly larger than the one employed in GLIMSE, with aperture corrections derived based on fluxes measured through Kron ellipses and the growth curves of the PSFs. Photometric uncertainties were estimated by placing random apertures in blank sky regions, accounting for correlated noise.

# 4. Results

## 4.1. Pop III Candidates Found in the Deepest NIRCam Data

We systematically applied our fiducial selection method (Section 2.5) to the NIRCam photometry catalogs from GLIMPSE, UNCOVER, CEERS, PRIMER, and JOF. We also performed careful visual inspections of the NIRCam images for all candidates to exclude potential artificial contaminants, such as residual diffraction spikes or extended flux from nearby bright objects. As a result, we identified one plausible candidate in GLIMPSE, GLIMPSE-16043, which satisfies all selection criteria. In Tables 3 and 4, we summarize the NIRCam photometry and the χ² values from different SED fitting methods for GLIMPSE-16043. We also identified a tentative candidate in JOF, JOF-21739, which meets the SED thresholds but falls very close to the borders of the criteria in the color–color diagrams, failing to satisfy all three. For completeness, we further investigate and discuss JOF-21739 in Appendix D. However, several suspicious aspects remain in this object, and thus, JOF-21739 should be regarded as a tentative candidate.

In Figure 6, we present the NIRCam colors of GLIMPSE-16043 (filled pentagon) and JOF-21739 (open pentagon) in the color–color diagrams. For reference, we also show other NIRCam sources with z<sub>phot</sub> = 5.0–7.5 in GLIMPSE (gray dots) and three sources that meet the criteria of the color-based criteria but do not satisfy the SED-based criteria (black open circle). GLIMPSE-16043 resides within the color space of Pop III galaxies and is well separated from the region occupied by metal-enriched galaxies. This clear separation is also confirmed in the other three color–color diagrams with F480M (see Appendix B), ensuring its robustness as the candidate.

<footnote n="35">³⁵ https://jwst-uncover.github.io/megascience</footnote>

<footnote n="36">³⁶ https://dawn-cph.github.io/dja/imaging/v7/</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Same NIRCam color–color diagrams as Figure 2, but with actual sources. The red filled pentagon represents a promising Pop III candidate, GLIMPSE-16043 at $z \simeq 6.5$ found in GLIMPSE, which meets all three color–color diagrams and the SED criteria. The open red pentagon denotes a tentative candidate, JOF-21739 at $z \simeq 6.2$, found in JOF (see Appendix D). The colored dots represent 450 NIRCam sources in the GLIMPSE field with photometric redshifts in the range $z_{phot} = 5.0-7.5$, where the color scale indicates their redshift. The open black circles present those that meet the three color–color diagrams, while they do not satisfy the SED criteria and thus do not pass our selection. The ~1% contamination rate observed in the color-based method alone (see top-right panel of Figure 3), suggests that a few potential contaminations may exist among the 450 NIRCam sources at $z_{phot} = 5.0-7.5$. The number of these color-satisfying but non-SED-consistent sources agree with the prediction from the contamination rate estimate, highlighting the importance of combining the color- and SED-based methods.</figcaption></figure>

<table><caption>**Table 4**
SED Fitting Results for GLIMPSE-16043 with Different Codes</caption><thead><tr><th>Fitting Type</th><th>Pop III (EAZY)</th><th>Galaxy (EAZY)</th><th>Galaxy (BEAGLE)</th><th>Galaxy (Prospector)</th><th>Galaxy (BAGPIPES)</th></tr></thead><tbody><tr><td>z<sub>phot</sub></td><td>6.50<sup>+0.03</sup><sub>-0.24</sub></td><td>6.18<sup>+0.15</sup><sub>-0.03</sub></td><td>6.28<sup>+0.12</sup><sub>-0.09</sub></td><td>6.41<sup>+0.08</sup><sub>-0.10</sub></td><td>6.26<sup>+0.15</sup><sub>-0.11</sub></td></tr><tr><td>χ<sup>2</sup></td><td>2.7</td><td>13.1</td><td>7.7</td><td>5.9</td><td>7.6</td></tr></tbody></table>

JOF-21739 meets one of the color–color diagrams and falls slightly outside the thresholds in the other two diagrams, while still satisfying the thresholds within the $1\sigma$ errors.

Figure 4 shows the location of GLIMPSE-16043 in the completeness and contamination rate distributions derived from the combined color-based and SED-based selection functions. GLIMPSE-16043 is detected with F200W = 30.0 mag (observed frame) and has a photometric redshift of $z_{phot,PopIII} = 6.50_{-0.24}^{+0.03}$. At this parameter space, the completeness is estimated at 50%, and the contamination rate is 0%, supporting the robustness of the selection. Since JOF-21739 does not meet all of the selection criteria of our fiducial method, we rerun the completeness and contamination rate tests optimized to the criteria satisfied by JOF-21739. We confirm that JOF-21739 also achieves high completeness (>80%) and zero contamination rate at its specific parameter space, while a modest contamination rate (~1%) remains in the nearby parameter space, making it a tentative Pop III candidate.

Figure 7 shows NIRCam cutouts of GLIMPSE-16043. The faintness in F410M is visually confirmed, supporting the presence of Pop III-like key SED features (Figure 1). GLIMPSE-16043 appears as a moderately magnified source with a lensing magnification of $\mu = 2.9_{-0.2}^{+0.1}$ with radial × tangential magnifications of $1.28 \times 2.26$, based on an updated mass model incorporating the GLIMPSE data (L. Furtak et al. 2025, in preparation). After lensing correction, the intrinsic absolute UV magnitude of GLIMPSE-16043 is $M_{UV} = -15.89_{-0.14}^{+0.12}$, indicating an ultrafaint source. The successful identification of candidates in these faintest UV regimes aligns with the predicted optimal depths for detecting galaxies dominated by Pop III stars (Section 2.6).

Figure 7 also presents the best-fit SEDs for GLIMPSE-16043 and metal-enriched galaxy templates using EAZY, the methodology introduced in Section 2.2. The photometric

redshift solutions are securely constrained by the Lyman break in the NIRCam+HST filters at ~0.9 µm and the enhanced flux observed in F444W+F480M, consistent with Hα emission at the inferred Lyman-break redshift, while possible low-z solutions are further discussed in Section 5.3. The HST +NIRCam SED is well fitted by a combination of the young (1–50 Myr) Pop III.2 models (Section 2.2) with $\chi^2 = 2.7$ and no dust attenuation, yielding a $\Delta\chi^2 = 10.4$ compared to the best-fit EAZY metal-enriched galaxy SED. This large $\Delta\chi^2$ value confirms that GLIMPSE-16043 differs significantly from typical metal-enriched galaxies.

To investigate whether the HST+NIRCam SED is truly challenging for metal-enriched galaxy models beyond EAZY, we further conducted SED fitting using the more flexible SED fitting tools of BEAGLE (J. Chevallard & S. Charlot 2016), BAGPIPES (A. C. Carnall et al. 2018), and Prospector (B. D. Johnson et al. 2021). We assumed parameter ranges typical for high-z star-forming galaxies (e.g., N. A. Reddy et al. 2023; M. Tang et al. 2023) with the redshift range $z = 0–8$. The detailed setup of the BEAGLE fits is as described in R. Endsley et al. (2025); Prospector fits are as described in R. P. Naidu et al. (2024); and BAGPIPES fits are as described in J. Chisholm et al. (2025, in preparation). Figure 8 summarizes the best-fit SEDs obtained from these tools. While the $\Delta\chi^2$ values decrease compared to EAZY, they generally remain >3–5 in all cases, reinforcing a preference for the Pop III model. It is worth mentioning that all best-fit SEDs reach the lowest gas-phase metallicity values in their grids ($Z_{gas}/Z_{\odot} \simeq 0.01$), suggesting that even if our candidates are not true Pop III galaxies, they likely represent extremely metal-poor galaxies. We further discuss alternative scenarios, including the extremely metal-poor galaxy case, in Section 5.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Best-fit SED of GLIMPSE-16043 derived with *EAZY* using the Pop III templates (blue curve) and the enriched galaxy templates (gray curve), together with the NIRCam and HST 3" × 3" cutouts in the top panels. The open red circles denote the NIRCam and HST/F814W photometry, and the red arrows indicate 2σ upper limits. The blue and gray squares present the expected photometry from the best-fit SEDs. The photometry favors the best-fit Pop III templates. The complete absence in F410M strongly supports the presence of an extremely strong Hα (EW ≈ 3000Å) and a Balmer jump that are consistent with Pop III SEDs but are challenging to reproduce along with the other emission-line excesses in metal-enriched galaxy SEDs.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Best-fit metal-enriched galaxy SEDs generated using *BAGPIPES* (green curve), *BEAGLE* (orange curve), and *Prospector* (magenta curve). The red symbols are the same photometric points as those shown in Figure 7. Despite the greater flexibility offered by these SED fitting codes compared to the Pop III template fitting with *EAZY*, the Δχ² values still exceed 3–5 in all cases. Additionally, all models reach the lowest Zgas values within their respective grids, Zgas/Z⊙ ≈ 0.01, where zero metallicity remains consistent within the uncertainties.</figcaption></figure>

## 4.2. Physical Properties

We summarize the physical properties of GLIMPSE-16043 in Table 5. Below, we briefly describe how these properties were derived.

<table><caption>**Table 5**
Physical Properties of GLIMPSE-16043</caption><thead><tr><th>Property</th><th>Measurement</th></tr></thead><tbody><tr><td>z<sub>phot</sub></td><td>6.50<sup>+0.03</sup><sub>-0.24</sub></td></tr><tr><td>μ</td><td>2.9<sup>+0.1</sup><sub>-0.2</sub></td></tr><tr><td>β</td><td>-2.34 ± 0.36</td></tr><tr><td>M<sub>UV</sub> [mag]</td><td>-15.89<sup>+0.12</sup><sub>-0.14</sub></td></tr><tr><td>M<sub>star</sub> [M<sub>☉</sub>]</td><td>≈10<sup>5</sup></td></tr><tr><td>r<sub>e</sub> [pc]</td><td>&lt;40</td></tr><tr><td>t<sub>age</sub> [Myr]</td><td>2.8</td></tr><tr><td>Hα EW [Å]</td><td>2810 ± 550</td></tr><tr><td>OIII/Hβ<sup>a</sup></td><td>&lt;0.44</td></tr><tr><td>12 + log(O/H)<sup>a</sup></td><td>&lt;6.4</td></tr></tbody></table>

**Note.** Physical parameters are derived using the best-fit Pop III templates (see Section 4.2 for details).

<sup>a</sup> The 1σ upper limit is presented, same as E. Vanzella et al. (2023b).

### 4.2.1. EW(Hα), Stellar Age, and UV Continuum Slope

The significant faintness in F410M observed in GLIMPSE-16043 strongly supports the presence of a Balmer jump and strong Hα emission. By determining the underlying continuum with the best-fit Pop III SEDs and utilizing the flux excess in F480M and F444W, we infer the rest-frame Hα equivalent width (EW) of 2810 ± 550 Å. These EW(Hα) values are comparable to or even higher than the maximal EW(Hα) values covered by the Yggdrasil models (E. Zackrisson et al. 2011), indicating that the majority of the emission is
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The point-source morphology of the Pop III candidate is consistent with arising from, e.g., extremely compact star clusters. NIRCam F150W 1" × 1" cutout for the observed, PSF, and PSF-subtracted residual image of GLIMPSE-16043. The compactness is consistent with the PSF, and we obtain an upper limit of 40 pc for the effective radius after the lens correction.</figcaption></figure>

dominated by very young stellar populations (<5 Myr; see, e.g., Figure 9 in J. A. A. Trussler et al. 2023) and/or the presence of very massive stars (>100 $M_\odot$) in these objects (D. Schaerer et al. 2025). From the best-fit Pop III SEDs, we obtain the mean $t_{age}^{37}$ of 2.8 Myr, which is consistent with the very young stellar ages implied by the high EW(H$\alpha$) above. Using the rest-frame UV filters F115W, F150W, and F200W, we also estimate a UV continuum slope of $\beta = -2.34 \pm 0.36$. This value is consistent with the nebular-dominated UV continuum shape predicted for Pop III SEDs in young stellar ages (~3 Myr) with no dust attenuation (see, e.g., Figure 7 in J. A. A. Trussler et al. 2023), well aligning with the above $t_{age}$ estimate and the observed high EW(H$\alpha$).

### 4.2.2. Stellar Mass

In the Yggdrasil models, the stellar component is uniformly set to $M_{star} = 10^6 M_\odot$. By scaling each template to the best-fit Pop III SED, we obtain $M_{star}$ values of $3.3 \times 10^5 M_\odot$ after the lens correction. We note, however, that the SED shapes of our Pop III candidates are dominated by the nebular continuum. This indicates that the $M_{star}$ estimates for our candidates strongly depend on the ionizing properties of galaxies, which are regulated by the IMF shape, $t_{age}$, and the nebular conditions such as electron density and ionization parameter. For example, if we simply assume the Pop III.1 IMF and $t_{age} = 1$ Myr, given the high EW(H$\alpha$) of ≈2800 Å, the $M_{star}$ estimates decrease to $0.7 \times 10^5 M_\odot$. While a detailed investigation of the ionizing properties including nebular conditions in Pop III galaxies is beyond the scope of this work, these results indicate that a simple order-of-magnitude estimate for $M_{star}$ would be ≈$10^5 M_\odot$ for GLIMPSE-16043, which is in excellent agreement with the baseline $M_{star}$ expected for newly forming Pop III galaxies at $z \simeq 6-7$ (Section 2.6).

### 4.2.3. Size

GLIMPSE-16043 exhibits extremely compact morphology. Using empirical PSFs in GLIMPSE generated in the same manner as A. Weibel et al. (2024), we perform Sérsic profile fitting with GALFIT (C. Y. Peng et al. 2010) on the high-resolution F150W image. Notably, the fitting does not converge, reaching the smallest size in the fitting grid, indicating consistency with a point source. Figure 9 presents

the NIRCam F150W 1" × 1" map around GLIMPSE-16043, the PSF in F150W, and the residual map after PSF subtraction. The morphology closely matches the PSF, supporting its compact nature. Typically, intrinsic sizes can be reliably constrained down to ~0.5 times the pixel scale (see Appendix C in M. Messa et al. 2022). Additionally, Y. Ono et al. (2023) conducted Monte Carlo simulations and found that the output/input size ratio for the smallest sources (effective radius $r_e = 1$ pix) with luminosities near the detection limit is approximately ~2 in NIRCam/SW filters (see Figure 4 and 5 in Y. Ono et al. 2023). By conservatively adopting the possible output/input size ratio of 3, we place an upper limit on the effective radius of GLIMPSE-16043 at $0.5 \times$ pixel scale (=0.02) × 3 ≃ 40 pc after applying the lensing correction.38 Compact star-forming clumps with $r_e < 10–100$ pc have been identified in recent JWST studies of strongly lensed galaxies at $z \gtrsim 6$ within similarly UV-faint and low-mass regimes (e.g., E. Vanzella et al. 2023b; A. Adamo et al. 2024; S. Fujimoto et al. 2024a; L. Mowla et al. 2024). The comparably small sizes of GLIMPSE-16043 may suggest it is a UV-faint, low-mass compact star-forming region or galaxy potentially dominated by Pop III stars.

### 4.2.4. Gas-phase Metallicity

In Figure 10, we present constraints on the [O III]/H$\beta$ ratio and the resulting gas-phase metallicity ($Z_{gas}$) of GLIMPSE-16043, derived from H$\alpha$ and [O III]+H$\beta$ fluxes inferred through flux excesses in the F480M (or F444W) and F356W filters, respectively. The best-fit Pop III SED model is used to estimate the underlying continuum and calculate the flux excesses. Although the F480M (or F444W) filter captures both H$\alpha$ and [N II] lines, the [N II]/H$\alpha$ ratio is expected to be ≲1% in the low-$Z_{gas}$ regime (e.g., Y. I. Izotov et al. 1997), making the contribution of [N II] to the flux in F480M or F444W negligible, if any. The excellent agreement of GLIMPSE-16043 with the Pop III templates in the rest-frame UV suggests nearly zero dust obscuration. Consequently, we estimate the H$\beta$ line strength by assuming an intrinsic H$\alpha$/H$\beta$ ratio of 2.74, with Case B recombination under conditions of an electron temperature of $2 \times 10^4$ K and an electron density of $10^4$ cm$^{-3}$ (D. E. Osterbrock 1989). From this analysis, we find the [O III] line flux to be consistent with zero and measure an H$\beta$ line flux of $(2.5 \pm 0.5) \times 10^{-19}$ erg s$^{-1}$ cm$^{-2}$. This yields a $1\sigma$ upper limit of 0.36, derived from the propagated uncertainties, which are primarily limited by the current NIRCam imaging depth.

The left panel of Figure 10 extrapolates photoionization modeling to the lowest $Z_{gas}$ regime along with a log(U) grid. Given the observed large EW(H$\alpha$) values in our candidates, we adopt log(U) ≳ −1.5 as the most plausible scenario and estimate upper limits of 12 + log(O/H) at ≲6.4 (≈0.005Z$_\odot$), based on the [O III]/H$\beta$ upper limits derived above. These constraints represent some of the lowest $Z_{gas}$ values observed to date at $z > 6$, comparable to LAP1 and LAP1-B, the most metal-poor clumps identified thus far, found within a strongly lensed arc at $z = 6.6$ (E. Vanzella et al. 2023a). A critical distinction is that LAP1 and LAP1-B lack continuum detections. This suggests that they may represent extremely metal-poor and low-mass (≈$10^3 M_\odot$) star clusters inside the

<footnote n="38">38 In this system, tangential magnification dominates and is used for lens correction rather than a circularized value.</footnote>

<footnote n="37">37 Luminosity-weighted logarithmic mean.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Constraints on the gas-phase metallicity $Z_{\text{gas}}$. Left panel: the red pentagon and the upper error bar represent the best-fit estimate and the $1\sigma$ upper limit on the [O III]/H$\beta$ line ratio, derived from the flux excesses in F480M (for H$\alpha$) and F356W (for [O III]+H$\beta$). In this calculation, the underlying continuum is evaluated using the best-fit Pop III template, and H$\alpha$/H$\beta$ is fixed at 2.74 (see Section 4.2.4). The dashed curves indicate photoionization model predictions extrapolated toward the low-metallicity regime from E. Vanzella et al. (2023a), covering an ionization parameter range of [-0.5: -2.0]. The red diagonal arrow shows how the upper limit of 12 + log(O/H) shifts according to that of [O III]/H$\beta$. The green symbols present the very [O III]-weak H$\alpha$ emitters at $z = 6.6$, LAP1 and LAP1-B (E. Vanzella et al. 2023b) found in a strongly lensed arc with no continuum detection, suggesting their possible origins in nebular emission arising from metal-free pockets within galaxies. Right panel: luminosity–metallicity relation. Given the very high EW(H$\alpha$) of ≈3000 Å, we assume log(U) ≥ -1.5 for GLIMPSE-16043 to derive the metallicity from the photoionization model shown in the left panel. The gray symbols denote the recent JWST measurements for lensed (I. Chemerynska et al. 2024) and field galaxies (K. Nakajima et al. 2023) at $z \simeq 4-10$, and the gray shaded region indicates the $1\sigma$ of its best-fit linear relation. LAP1 and LAP1-B fall on the best-fit relation. In contrast, GLIMPSE-16043 shows significant UV enhancement relative to the best-fit relation, which indicates a unique origin, e.g., due to a top-heavy IMF in a Pop III galaxy or an extremely metal-poor, low-mass AGN. Under the condition that the dust cooling is minimal, [O/H] ≃ -3 is a critical value for the Pop III to Pop II transition via metal cooling (e.g., V. Bromm & A. Loeb 2003), indicating that $[Z_{\text{gas}}/Z_{\odot}] < -3$ is the Pop III regime (cyan shade).</figcaption></figure>

metal-enriched galaxy or line-emission from pristine nebular gas clouds within the galaxy, rather than star clusters. In contrast, our candidates potentially represent the lowest metallicity observed as a single, isolated galaxy at $z > 6$.

The right panel of Figure 10 compares the luminosity–metallicity relation of our candidates with recent JWST measurements of $z \simeq 4-10$ galaxies (e.g., K. Nakajima et al. 2023; I. Chemerynska et al. 2024) and a linear fit to these observations (black shaded region). While LAP1 and LAP1-B fall along the best-fit relation, GLIMPSE-16043 and JOF-21739 exhibit significant deviations, appearing more than 100 times brighter in the rest-UV than LAP1 and LAP1-B at comparable metallicities. This deviation may indicate mechanisms that enhance UV luminosity at a given stellar mass, such as a top-heavy IMF, which is consistent with theoretical predictions for Pop III galaxies (e.g., J. Tumlinson 2006; M. A. Fardal et al. 2007; R. Davé 2008). Alternatively, the deviation could be attributed to significantly lower $Z_{\text{gas}}$ at a given $M_{\text{UV}}$, potentially resulting from a recent pristine gas inflow that dilutes the metal abundance and decreases the 12 + log(O/H) value. We further explore possible mechanisms behind this deviation from the typical luminosity–metallicity relation in Sections 4.5 and 5.

Although the definition of Pop III stars is strictly zero metallicity, it is useful to consider a critical metallicity threshold ($Z_{\text{crit}}$), above which “normal” (low-mass dominated) star formation can occur (i.e., the transition from Pop III to Pop II). This helps us understand how close our candidates may be to this boundary. Theoretically, metal-line cooling can set $Z_{\text{crit}}$ at about [O/H] ~ -3 (V. Bromm & A. Loeb 2003), whereas dust

cooling could drive the threshold even lower, to [O/H] ~ -5 or -6 (R. Schneider et al. 2006). However, PISNe from very massive primordial stars can produce extremely energetic explosions, often resulting in significant destruction of newly formed dust grains (e.g., T. Nozawa et al. 2010). Furthermore, once the first supernovae occur, metal enrichment within the star-forming regions proceeds rapidly, so that [O/H] ≈ -3 can be reached on short timescales (≈10 Myr; e.g., T. Karlsson et al. 2008; U. Maio et al. 2010). Given uncertainties in dust survival and the rapid pace of chemical enrichment, we adopt [O/H] = -3 as a practical reference for $Z_{\text{crit}}$. In Figure 10, the region below $Z_{\text{crit}}$ is shaded in cyan to indicate the Pop III regime. Our candidate, with [O/H] < -2.3 falls within a factor of a few from $Z_{\text{crit}}$. Future deep observations may push their upper limits farther into the Pop III regime.

### 4.3. Pop III UVLF at $z \simeq 6-7$

The presence of late Pop III star formation and Pop III galaxies down to $z \simeq 6$ largely depends on several factors, including feedback mechanisms from early galaxies, the efficiency of subsequent chemical mixing, environmental conditions in star-forming regions, and the abundance of the metal-free low-mass halos remaining at $z \simeq 6$. Conversely, the volume density of Pop III galaxies at $z \simeq 6$ provides critical constraints for simulations of Pop III formation and evolution in the early Universe (e.g., A. Pallottini et al. 2014a; R. Sarmento et al. 2018; B. Liu & V. Bromm 2020; E. Visbal et al. 2020; A. Vikaeus et al. 2022; A. Venditti et al. 2023, 2024b). Using the systematic search for Pop III galaxies with NIRCam and the well-defined selection functions
developed in this study (Section 2), we provide the first observational constraint on the Pop III UV luminosity function (UVLF). Although spectroscopic confirmation is essential to definitively determine whether these candidates are true Pop III galaxies, we derive the Pop III UVLF at $z \simeq 6-7$ using our photometric candidates, following the standard approach commonly employed in UVLF studies of high-redshift galaxies. The dedicated survey volume, including completeness corrections, also allows us to set meaningful upper limits even if future spectroscopy reveals these candidates to be enriched sources.

First, we evaluate the survey volume as follows. Based on the completeness function shown in Figure 3, we focus on the parameter space satisfying the completeness of $\ge 50\%$, defining our effective search range as $z = 5.6-6.6$ and observed magnitudes in F200W of $\le 30.3$ mag. Note that the completeness in Figure 3 is estimated using GLIMPSE data. When considering completeness in other survey fields, we scale the distribution according to the relative data depth in each survey compared to GLIMPSE, while the completeness distribution remains largely consistent across redshifts due to the similar NIRCam filter sets. For GLIMPSE, we calculate the effective survey area using the NIRCam field of view of 9.4 arcmin² and the updated mass model incorporating the GLIMPSE data (L. Furtak et al. 2025, in preparation; I. Chermynska et al. 2025, in preparation). For UNCOVER, the total NIRCam mosaic area is 45 arcmin² (R. Bezanson et al. 2024), and we also calculate the effective survey area using the publicly available v2.0 mass model<fnref n="39" /> (L. J. Furtak et al. 2023; S. H. Price et al. 2024; K. A. Suess et al. 2024; J. R. Weaver et al. 2024). For CEERS and JOF, we adopt areas of 88.1 arcmin² (S. L. Finkelstein et al. 2024) and 9.4 arcmin² (B. Robertson et al. 2024), respectively. For PRIMER, we follow the areas defined by the depths in C. T. Donnan et al. (2024) and adopt values of 47.8 arcmin², 84.4 arcmin², and 170.7 arcmin² for PRIMER-COSMOS Deep, PRIMER-COSMOS Wide, and PRIMER-UDS, respectively. We then bin the effective survey areas according to the absolute UV magnitude ($M_{UV}$) bins, resulting in a total survey area of approximately 500 arcmin² in the brightest bin. We convert the survey area into a comoving survey volume for the $z = 5.6-6.6$ slice and apply completeness corrections to our individual detections and upper limits for the final volume density estimates. Uncertainties and upper limits are calculated using Poisson errors following the values presented in N. Gehrels (1986). Note that the observed magnitude and intrinsic $M_{UV}$ are not directly correlated due to varying magnifications in lensing fields, and we may calculate the completeness for each magnification at a given observed magnitude. However, our completeness simulations assume point sources (i.e., unaffected by lensing distortions), making the completeness depend solely on the magnitude in the observed frame. We thus adopt an average completeness of 75% based on the effective survey parameter space defined above in the completeness correction for the upper-limit estimates in lensing fields, while we use the observed magnitude and redshift estimate and infer the completeness for our individual candidates.

In Figure 11, we show the inferred Pop III UVLF at $z \simeq 6-7$, derived from our robust candidate GLIMPSE-16043 (red filled

pentagon). For completeness, we also plot the tentative candidate JOF-21739 (red open pentagon). Colored curves and shaded regions represent simulation predictions from the literature (B. Liu & V. Bromm 2020; E. Visbal et al. 2020; R. Sarmento & E. Scannapieco 2022; A. Venditti et al. 2023, 2024b), updated here to match the UVLF format (see Appendix E). Our Pop III UVLF estimates, including the upper limits, are listed in Table 6, and the theoretical framework for each Pop III UVLF simulation prediction is summarized in Appendix E.

We find that GLIMPSE-16043 aligns remarkably well with several simulation models, despite considerable differences in the assumptions for both galaxy formation and evolution and their implementations of the UV emission (see details in Appendix E). This lends additional support to the scenario that GLIMPSE-16043 is indeed a Pop III galaxy, independently validated from the perspective of volume density. On the other hand, most of these simulations exhibit a sharp cutoff beyond $M_{UV} \simeq -16$, corresponding to a mass cutoff around $M_{star} \approx 10^5 M_{\odot}$ due to the limited survey volumes in the simulations and the specific conditions required for Pop III star formation at $z \sim 6-7$ (Section 2.6). Note that this does not immediately rule out the presence of JOF-21739, given the significant uncertainties in the theoretical upper limit of the Pop III galaxy mass (e.g., H. Yajima & S. Khochfar 2017) and in the mass-to-light ratio of Pop III galaxies. If physical mechanisms enhance UV luminosity at a given stellar mass, the current shapes of the simulated Pop III UVLFs may extend farther into the bright end. Possible mechanisms include rapidly rotating stars undergoing chemically homogeneous evolution, which could increase the brightness by up to $\sim 2$ mag (B. Liu et al. 2025), thus making the current Pop III UVLF predictions consistent also with JOF-21739.

### 4.4. Pop III SFRD at $z \simeq 6-7$

We also explore possible constraints on the cosmic Pop III star formation rate density (SFRD), $\rho_{*,PopIII}$, based on our UVLF results. First, for the lower limit, we focus on the single data point of GLIMPSE-16043 to obtain a conservative estimate. We use its lower-limit volume density to derive the corresponding UV luminosity density. D. Schaerer et al. (2025) examined the variation of the conversion factor between UV luminosity and SFR, $\kappa_{UV}$,<fnref n="40" /> for different IMFs and stellar populations (including Pop III) as a function of $t_{age}$ (see Figure 7 in D. Schaerer et al. 2025). Based on the best-fit mean $t_{age} = 2.8$ Myr and a moderately top-heavy IMF in GLIMPSE-16043, we adopt $\kappa_{UV} = 0.5 \times 10^{-28}$ and obtain $\rho_{*,PopIII} = 1.5 \times 10^{-6} [M_{\odot} \text{yr}^{-1} \text{cMpc}^{-3}]$.

Second, for the upper limit, we use the best-fit UVLF shape at $z = 6$ derived in R. J. Bouwens et al. (2022), scaling it to match the upper limit of the UVLF data point of GLIMPSE-16043. We then integrate the rescaled $z = 6$ UVLF down to $M_{UV} = -12$, which reflects the possible $M_{UV}$ range of Pop III galaxies down to $\approx 10^4 M_{\odot}$. Using the same $\kappa_{UV}$ value, we derive $\rho_{*,PopIII} = 1.5 \times 10^{-4} [M_{\odot} \text{yr}^{-1} \text{cMpc}^{-3}]$. As discussed in Section 4.3, there is likely a cutoff in the bright end of the Pop III UVLF based on simulation predictions. Nevertheless, using the same UVLF shape as the general galaxy population provides a conservative upper-limit estimate. We summarize our possible Pop III SFRD constraints in Table 7. Note that we

<footnote n="39">³⁹ https://jwst-uncover.github.io/DR4.html#LensingMaps</footnote>

<footnote n="40">⁴⁰ Defined as per SFR(UV) [$M_{\odot} \text{yr}^{-1}$] = $\kappa_{UV} \times L_{UV}[\text{erg s}^{-1} \text{Hz}^{-1}]$.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Pop III UV luminosity function (UVLF) at $z \simeq 6-7$. The red symbols represent the volume density estimates of our Pop III candidates identified through a systematic search across a total area of $\simeq 500 \text{ arcmin}^2$ from the JWST legacy surveys of GLIMPSE, UNCOVER, CEERS, PRIMER, and JOF. Filled and open pentagons correspond to GLIMPSE-16043 and JOF-21739, respectively, with error bars reflecting Poisson uncertainties. Upper limits are calculated based on Poisson uncertainties at the single-sided confidence level of 84.13% (N. Gehrels 1986). The black curve shows the UVLF of general $z \simeq 6-7$ galaxies for comparison (R. J. Bouwens et al. 2022), while the colored curves and shades represent theoretical predictions (see details in Appendix E). GLIMPSE-16043 falls well within the range of theoretical predictions in the ultrafaint regime ($M_{UV} \simeq -16$), whereas JOF-21739, if real, may require additional physical mechanisms, such as rapidly rotating stars (B. Liu et al. 2025), to explain its boosted $M_{UV}$ for some models.</figcaption></figure>

<table><caption>**Table 6**
Constraints on Pop III UVLF at $z \sim 6-7$</caption><thead><tr><th>M<sub>UV</sub><br>(AB mag)</th><th>Φ<br>(10<sup>-5</sup> Mpc<sup>-3</sup> dex<sup>-1</sup>)</th></tr><tr><th>(1)</th><th>(2)</th></tr></thead><tbody><tr><td>-19.5</td><td>&lt;0.16</td></tr><tr><td>-18.5</td><td>&lt;0.35</td></tr><tr><td>-17.5</td><td>(0.64<sup>+1.46</sup><sub>-0.53</sub>)<sup>a</sup></td></tr><tr><td>-16.5</td><td>&lt;5.0</td></tr><tr><td>-15.5</td><td>18.0<sup>+41.4</sup><sub>-14.0</sub></td></tr><tr><td>-14.5</td><td>&lt;226.5</td></tr></tbody></table>

**Notes.** (1): Absolute UV magnitude. (2): Observational constraints obtained from GLIMPSE and other JWST legacy deep field data of UNCOVER (R. Bezanson et al. 2024), PRIMER-UDS/COSMOS (C. T. Donnan et al. 2024), and CEERS (S. L. Finkelstein et al. 2024), resulting in a total survey area of $\simeq 500 \text{ arcmin}^2$ (0.14 deg$^2$). Errors and upper limits are $1\sigma$, evaluated with the Poisson uncertainty (N. Gehrels 1986).

<sup>a</sup> This is estimated from the tentative candidate of JOF-21739. If we place an upper limit instead, the $1\sigma$ upper limit is estimated to be $1.2 \times 10^{-5} \text{ Mpc}^{-3} \text{ dex}^{-1}$.

focus on GLIMPSE-16043 for our Pop III SFRD estimates above, because of the less robust nature of JOF-21739. However, the contribution of JOF-21739 still falls well within our lower and upper limits, given our conservative approach.

In Figure 12, we show our possible constraints on the cosmic Pop III SFRD (red line and shading). For reference, we also show the total SFRD (black curve) by integrating the best-fit redshift evolution of the general galaxy UVLFs from R. J. Bouwens et al. (2022) down to $M_{UV} = -12$, the same

<table><caption>**Table 7**
Constraints on Pop III SFRD at $z \sim 6-7$</caption><thead><tr><th>z<br>(1)</th><th>log(ρ<sub>*</sub>,PopIII)<br>(M<sub>☉</sub> yr<sup>-1</sup> Mpc<sup>-3</sup>)<br>(2)</th></tr></thead><tbody><tr><td>5.6–6.6</td><td>[-5.81; -3.82]</td></tr></tbody></table>

**Note.** (1): Redshift range (see Section 4.3). (2): Our lower- and upper-limit estimates by accounting for the single galaxy contribution of GLIMPSE-16043 and integrating the $z = 6$ UVLF for general galaxies (R. J. Bouwens et al. 2022) rescaled down to match the data point of GLIMPSE-16043 (see Section 4.4).

limit used for the Pop III calculation. For comparison, we include simulation predictions from the literature (A. Pallottini et al. 2014b; U. Maio et al. 2016; H. Xu et al. 2016; J. Jaacks et al. 2019; B. Liu & V. Bromm 2020; D. Skinner & J. H. Wise 2020; E. Visbal et al. 2020; J. B. Muñoz et al. 2022; R. Sarmento & E. Scannapieco 2022; S. Hegde & S. R. Furlanetto 2023; A. Venditti et al. 2023). Our Pop III SFRD constraint at $z \simeq 6-7$ agrees well with the broad range of simulation predictions. Similar to our UVLF results, this SFRD constraint further supports the robust identification of GLIMPSE-16043 as a Pop III candidate. Compared to the total SFRD, our results suggest that Pop III galaxies may account for ~0.01%–1% of the total cosmic SFRD at $z \simeq 6-7$.

We caution that contributions from metal-free pockets/satellites, (i.e., Pop II+Pop III mixed objects; candidates—e.g., E. Vanzella et al. 2023b; X. Wang et al. 2024) are missed in our current Pop III SFRD estimates, due to the strict requirement for the absence of the [O III] line in our selection
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Cosmic Pop III star formation rate density (SFRD) at $z \simeq 6-7$. The red shading and line indicate our conservative lower- and upper-limit estimates, respectively, derived from our Pop III UVLF results (see the text). Colored curves represent theoretical predictions for the Pop III SFRD from the literature (A. Pallottini et al. 2014b; U. Maio et al. 2016; H. Xu et al. 2016; J. Jaacks et al. 2019; B. Liu & V. Bromm 2020; D. Skinner & J. H. Wise 2020; E. Visbal et al. 2020; J. B. Muñoz et al. 2022; R. Sarmento & E. Scannapieco 2022; S. Hegde & S. R. Furlanetto 2023; A. Venditti et al. 2023). The black curve shows the total SFRD inferred from the best-fit redshift evolution of the general galaxy UVLFs in R. J. Bouwens et al. (2022), integrated down to $M_{UV} = -12$ to match a typical mass limit ($\sim 10^4 M_{\odot}$) in the simulations above. Our observational constraints fall well within the theoretical range and suggest that the Pop III galaxies may contribute ~0.01%–1% of the total cosmic SFRD at $z \simeq 6-7$.</figcaption></figure>

method (Section 2). We note, however, that we also applied the same Pop III selection technique to the clump-based NIRCam photometry for individual clumps in 133 strongly lensed, resolved galaxies at $z \simeq 5-8$ in UNCOVER and GLIMPSE (A. Claeyssens et al. 2025, 2025, in preparation), which resulted in null detection. Hence, the contribution from the PopII+PopIII mixed objects, even if any, would be unlikely to change our current SFRD constraints significantly. Concerted efforts will be required to derive a definitive conclusion about this fraction—e.g., the development of future JWST lensing cluster surveys.

We note that our selection is most sensitive to young Pop III populations ($\lesssim 20$ Myr; see Figure 2), but this age bias is not expected to significantly affect the derived SFRD constraints, as top-heavy IMFs are generally in pristine environments, resulting in short-lived massive stars that dominate the Pop III contribution (see Appendix E). However, if long-lived Pop III stars are more prevalent than currently predicted, both our observational constraints and the corresponding simulation predictions would increase accordingly, due to completeness corrections. This possibility highlights the need for future observations to better characterize the full age distribution of Pop III populations.

## 4.5. Cosmological Context of GLIMPSE-16043

The proposed selection technique in Section 2 and the Pop III candidate GLIMPSE-16043 provide a unique opportunity to study the pristine environments of early galaxies. Beyond the low implied $[O\,III]/H\beta$ ratio, GLIMPSE-16043 has an extreme $H\alpha$ equivalent width, blue continuum slope, and strong Balmer Jump signatures, suggestive of extremely

young stellar populations. Further, the relative faint $M_{UV}$ suggests that GLIMPSE-16043 resides in a low-mass halo. As such, GLIMPSE-16043 is significantly off the typical luminosity–metallicity relation. What do the observed properties of GLIMPSE-16043 suggest about the transformation of the pristine Universe into a near-universally metal-enriched Universe (e.g., T. Karlsson et al. 2013)? More concretely, how do galaxies enrich from the nearly pristine Pop III metallicity and onto the well-established luminosity–metallicity relation? As galaxies enrich, they must exceed the critical metallicity, $Z_{crit}$, where star formation transitions from the high-mass dominated Pop III mode to the low-mass dominated Pop II mode,<fnref n="41" /> but this exact transition has so far been observationally under-constrained.

In broad terms, there are two overall pathways for a star-forming system to achieve this transition: enriching above $Z_{crit}$ in one step from the first round of SNe, or from a multistep process. In general, the single enrichment scenario occurs very quickly (5–10 Myr), consistent with the extreme equivalent widths observed in GLIMPSE-16043. In the multiple episode scenario, enrichment from Pop III stars forming at $z \gtrsim 15$ in mini-halos establishes a seed metallicity that fuels subsequent star formation. The level of this pre-enrichment may still be below $Z_{crit}$, depending on the Pop III SN yields from a small number of massive stars and the mixing and dilution of their ejecta in the surrounding medium (e.g., M. Jeon et al. 2021). Since the pre-enriched gas is still below $Z_{crit}$, star formation requires subsequent accretion of gas and proceeds quite slowly

<footnote n="41">⁴¹ The precise value of $Z_{crit} \sim 10^{-6}-10^{-3} Z_{\odot}$ is under debate, involving dust- and fine-structure line cooling channels, with a possible dependence on redshift and local environment (e.g., R. Schneider et al. 2012).</footnote>
follow-up of the surrounding environment will also be crucial
once the Pop III nature of GLIMPSE-16043 is confirmed.

**5. Alternatives to the Pop III Scenario**

Although our Pop III candidates satisfy our stringent selection criteria and show excellent agreement in the volume density with the simulation predictions, which independently strengthen its robustness, it is important to explore all plausible alternative explanations given that we are searching for a very rare population. Below, we discuss interpretations for our candidates other than Pop III.

**5.1. Extremely Metal-poor Galaxy?**

The first plausible alternative scenario is an extremely metal-poor galaxy with $Z_{gas}/Z_{\odot}$ meeting the current upper limit of our candidates, but not entirely zero metallicity. This possibility naturally arises due to the limitation from the current photometric constraints on the [O III]/H$\beta$ ratio derived from the F356W excess. To draw a definitive conclusion, deep follow-up spectroscopy is essential to place more stringent upper limits on the complete absence of [O III] and/or to confirm the expected strong He II $\lambda$1640 and $\lambda$4686 equivalent width (J. Tumlinson et al. 2001; D. Schaerer 2002; K. Nakajima & R. Maiolino 2022).

Even if future spectroscopy reveals that our candidates are extremely metal-poor galaxies rather than Pop III, they still offer exciting new insights. First, the current upper limit ($Z_{gas}/Z_{\odot} < 0.005$) already reaches the most metal-poor regime explored by JWST at $z \gtrsim 6$ (Figure 10). This constraint is comparable to the values observed in the extremely [O III]-weak clumps in the strongly lensed arc at $z = 6.6$, known as LAP-1 and LAP-1B (E. Vanzella et al. 2023b). However, the rest-frame UV-optical continuum is undetected in LAP-1 and LAP-1B, leaving open the possibility that these clumps represent nebular emission from metal-poor or metal-free pockets of gas within galaxies rather than star clusters. In contrast, our candidates may serve as the first unique laboratory to study extremely metal-poor star clusters at $z > 6$. Second, the deviation in the luminosity-metallicity relation observed in our targets suggests either pristine gas inflow or enhanced UV flux. In the former case, directly detecting gas inflow has been challenging, particularly in ultrafaint, low-mass early galaxies with $M_{star} \approx 10^{5-6} M_{\odot}$. Our candidates provide key insights into the gas fueling of abundant low-mass nascent galaxies in the early Universe. In the latter case, the enhanced UV flux implies a high luminosity-to-mass ratio, offering compelling evidence of a top-heavy IMF in an extremely metal-poor galaxy at $z > 6$. Therefore, the extreme metallicity regime and the deviation in the luminosity-metallicity relation observed in our candidates represent in any case fairly unique features, even under the extremely metal-poor galaxy scenario.

**5.2. Extremely Metal-poor, Faint AGN?**

In the luminosity-metallicity relation (Figure 10), we find a significant deviation of our candidates from the typical relation, which could be attributed to either pristine gas inflow (i.e., reflecting an offset in $Z_{gas}$) or enhanced UV flux (i.e., reflecting an offset in $M_{UV}$). If the latter is true, possible mechanisms include a top-heavy IMF or AGN activity. While young star clusters may also exhibit compact morphologies

($\lesssim 100$ pc; e.g., A. Adamo et al. 2024; S. Fujimoto et al. 2024a; L. Mowla et al. 2024), the point-source morphology of our candidates is also consistent with the AGN scenario. Recent JWST studies have reported an abundant population of faint broad-line AGNs (e.g., D. D. Kocevski et al. 2023; J. Matthee et al. 2024; J. E. Greene et al. 2024; R. P. Naidu et al. 2024), with an estimated abundance of $\simeq 5\%–10\%$ among general galaxies at $z = 4–9$ down to $M_{UV} = -18$ (e.g., Y. Harikane et al. 2023b; S. Fujimoto et al. 2024c; R. Maiolino et al. 2024a; L. Napolitano et al. 2025). Interestingly, the volume density of our Pop III candidates is almost $\sim 1\%$ of the general galaxy population at $z \simeq 6–7$ (Figure 11), consistent with the abundance of faint AGNs, assuming their high abundance extends to fainter UV magnitudes. Taken together, these results suggest that the AGN scenario is worth further exploration as an alternative explanation. In particular, the UV luminosity of GLIMPSE-16043 ($M_{UV} = -15.9$) lies in an ultrafaint regime, even compared to the recent JWST studies of faint AGNs ($M_{UV} < -18$). Assuming a typical Eddington ratio of 0.1–1.0 observed in these JWST-discovered faint AGNs, this UV luminosity corresponds to a black hole (BH) mass of $M_{BH} \simeq 10^{4-5} M_{\odot}$. This mass range is consistent with seed BHs (e.g., M. Volonteri 2010; K. Inayoshi et al. 2020) and intermediate-mass BHs (e.g., J. E. Greene et al. 2020; P. Narajan 2021), which are thought to play a key role in the origins of supermassive black holes from the early Universe to the present day.

To explore this possibility, it is crucial to note that the extremely metal-poor conditions are still essential in the material surrounding the AGN, as evidenced by the absence of the [O III] line, which is typically very bright in standard SEDs of high-$z$ quasars and AGNs. We thus adopt several SED templates specifically constructed for a seed BH embedded in metal-poor environments. The first set is the model presented in K. Inayoshi et al. (2022b), assuming a seed BH mass of $M_{BH} = 10^6 M_{\odot}$ embedded in a metal-poor ($Z_{gas}/Z_{\odot} = 0.01$) galaxy. This SED includes three primary components: (1) radiation from the unresolved nuclear accretion disk surrounding the BH, modeled using a broken power-law spectrum (e.g., E. Lusso et al. 2015), (2) nebular emission lines and continuum reprocessed by the irradiated gas in the surrounding nebular region, calculated via CLOUDY simulations, and (3) radiation from the dense accretion disk (resolved in radiation-hydrodynamical simulations within 0.1–100 pc). The gas and disk boundary is defined where the electron fraction $x_e = 0.9$, marking the transition between the dense accretion disk and the surrounding nebular gas. The viewing angle is set to 60°.

The second set is the direct collapse BH (DCBH) model taken from K. Nakajima & R. Maiolino (2022). The model incorporates the following primary components: (1) an accretion disk emitting a thermal Big Bump component, parameterized by blackbody temperatures $T_{bb} = 5 \times 10^4$, $1 \times 10^5$, and $2 \times 10^5$ K, (2) a high-energy power-law continuum, described as $f_{\nu} \propto \nu^{\alpha}$, with spectral slopes $\alpha = -1.2, -1.6$, and $-2.0$, extending into the X-ray regime, and (3) nebular emission by using CLOUDY to simulate gas photoionized by the BH radiation with the plane-parallel geometry, $n_e = 10^3 \text{ cm}^{-3}$, $Z_{gas}/Z_{\odot} = 0-2$, $\log U = [-0.5: -3.0]$. The model does not include dust in the ionized regions, ensuring that emission lines are unaffected by dust absorption or depletion effects, which should not matter when considering extremely metal-poor or metal-free cases. For this study, we
Indeed, in Figure 14, we present color–color diagrams like
those used in our Pop III selection (Figure 2) for a compilation
of 2365 galaxies observed by JWST in the NIRSpec/prism
mode. Spectra of these sources were homogeneously reduced
with the msaexp software (G. Brammer 2022) and distributed
in the v3.0 release of the DAWN JWST archive (A. de Graaff
et al. 2025; K. E. Heintz et al. 2025).<sup>42</sup> We select galaxies with
the highest quality grade of 3 at z = 0–1.7, shift them to
z = 1.5, and synthesize NIRCam colors directly from the prism
spectra. Note that the color–color diagrams shown in Figure 2
are focused on LW bands, so we are able to use the lower-
redshift sources by ignoring the lack of prism coverage at bluer
wavelengths $\lesssim$0.7–0.8 µm. From this exercise, it is clear that
no source in the public JWST spectroscopic archive has colors
that match the color location of GLIMPSE-16043.

To summarize—the low-z scenario is still strongly dis-
favored by the forced low-z SED analysis and the NIRCam
color–color comparison with real DJA galaxies shifted at
z = 1.5, while the existence of a z ≃ 1.5 galaxy at r ∼ 4''6
warrants some caution. Importantly, these objects would have
to be extraordinary sources in their own right, the likes of
which JWST is yet to observe. Spectroscopic follow-up is
therefore imperative and the only way forward to confirm the
nature of these Pop III candidates.

5.4. Two Different Objects Mimicking Colors?

The chance projection of two objects at different redshifts
may alter the observed colors. Although our candidates exhibit
very compact morphologies, it is still useful to estimate the
probability of such a chance projection. In GLIMPSE, a total
of 59,650 sources is detected in the catalog, corresponding to a
surface density of NIRCam objects of 6345 arcmin⁻². Given
the point-source morphology of our candidates, we define a
projection radius of 0''01. Following the calculation presented
in A. J. B. Downes et al. (1986), we estimate the probability of
a chance projection to be ~0.06%. This scenario would also
require specific combined colors that satisfy all of the selection
criteria, including the color–color diagrams and SED thresh-
olds. As a result, the actual probability is likely reduced by
several orders of magnitude. We therefore conclude that this
scenario is highly unlikely.

5.5. What Is the Most Plausible Alternative Solution?

We briefly summarize the discussions in Sections 5.1–5.4. The extremely metal-poor galaxy scenario (Z<sub>gas</sub>/Z<sub>☉</sub> < 0.005) remains plausible, as it is consistent with the current photometric upper limit on the [O III]/Hβ ratio inferred from the F356W excess (Section 5.1). While the SED features observed in our Pop III candidates are difficult to reconcile with typical metal-poor AGN models, we cannot exclude the possibility of a metal-poor AGN due to the large uncertainties surrounding seed BH environments in early galaxies (Section 5.2). Contamination by low-z sources is strongly disfavored by the observed SED features, although we cannot rule out rare low-z objects that might exhibit similar SED properties yet remain undiscovered (Section 5.3). Finally, the chance projection of two distinct sources that collectively produce the observed peculiar colors is highly unlikely (Section 5.4).

Importantly, however, even under these alternate scenarios,
the faintness, the extremely low [O III]/Hβ ratio, and the
deviation in the M<sub>UV</sub>-Z<sub>gas</sub> relation observed in our candidates
suggest the identification of either metal-poor seed or
intermediate-mass black holes (≈10<sup>4</sup>-5 M<sub>☉</sub>), a recent pristine
gas inflow, or a top-heavy IMF formed in ultrafaint, low-mass
early galaxies with M<sub>star</sub> ≈ 10<sup>5</sup>-6 M<sub>☉</sub>. The low-z interlopers
would also be extreme starbursts in ultrafaint galaxies.
Regardless of their true origins—Pop III or otherwise—our
efficient NIRCam selection method opens a new discovery
space in the distant Universe.

**6. Future Prospects**

We have developed an efficient NIRCam-based method for
the Pop III galaxy search (Section 2), identified promising
candidates in the ~500 arcmin² area of publicly available deep
JWST legacy NIRCam data (Sections 3 and 4), and discussed
possible alternative explanations (Section 5). Here, we briefly
discuss future prospects based on these findings and the
current limitations.

First, deep follow-up spectroscopy is essential to confirm
the Pop III nature of the photometrically selected candidates.
Although our candidates show excellent agreement with
Pop III galaxy SED models, the presence of very weak [O III]
remains consistent with the uncertainty in the F356W
photometry, which naturally translates to the current upper
limits on [O III]/Hβ and subsequent Z<sub>gas</sub> estimates. This
limitation will be the same for all photometric candidates
identified in future surveys. To confirm the Pop III origin, an
alternative approach is the detection of strong He IIλ1640 and/
or He IIλ4686 EWs (e.g., V. Bromm et al. 2001; J. Tumlinson
et al. 2001; K. Nakajima et al. 2022; D. Schaerer 2002;
X. Wang et al. 2024). Thus, deep spectroscopy to place more
stringent upper limits on [O III]/Hβ and/or detect strong He II
is a logical next step. We recommend using NIRSpec medium-
or high-resolution gratings rather than the prism mode for such
follow-up observations. As discussed in Section 5, ruling out
the metal-poor AGN scenario is challenging. The use of
NIRSpec medium- or high-resolution gratings will enable the
detection of broad Balmer lines (e.g., Hα, Hβ) and confirm the
complete absence of neighboring [O III] lines and/or the large
EW of the neighboring He IIλ4686.

While the presence of a luminous z ~ 1.5 galaxy at r ~ 4''6
could be potentially concerning (Section 5.3), we must
remember that if such environmental considerations were
followed as an iron-clad rule, the most distant galaxy presently
known at z<sub>spec</sub> = 14.18, GS-z14-0 (S. Carniani et al.
2025, 2024; S. Schouws et al. 2025), would never have been
selected for follow-up observations. In particular, GS-z14-0's
closest neighbor, lying a mere <1''0 away, is a luminous
z = 3.47 galaxy with a redshift lining up exactly with the low-z
solution for GS-z14-0, with its Balmer break occurring at the
same wavelength as the Lyman break of the distant source
(S. Carniani et al. 2024). In the pursuit of extraordinary
sources, such as Pop III galaxy candidates, especially in the
early years of JWST, it is perhaps best to obtain follow-up
spectra liberally.

Second, observing a large number of lensing clusters is
critical for discovering (lensed) Pop III galaxy candidates and
enabling time-efficient follow-up spectroscopy. Assuming
M<sub>star</sub> ≈ 10<sup>5</sup> M<sub>☉</sub> as a typical stellar mass of the first galaxies
(Section 2.6), Figure 5 shows that a depth even fainter than

<footnote n="42"><sup>42</sup> See the Acknowledgments for a list of programs that these spectra were
collected as a part of.</footnote>
~30.5 mag is required, exceeding the deepest NIRCam survey depths achieved with more than 100 hr of observation time. This highlights the difficulty of discovering such objects in blank-field surveys. In this context, gravitational lensing is a promising way forward. The GLIMPSE program has demonstrated this potential by successfully identifying GLIMPSE-16043 within the survey volume of a single cluster. While observing a large number of lensing clusters comes at a cost of shallower depth per cluster, the depth is compensated by collecting numerous high-magnification patches. In fact, this “wide and shallow” strategy is the most time-efficient way to build up the deepest layer of imaging, as demonstrated analytically and by recent large lensing cluster observations (S. Fujimoto et al. 2024b; L. Vujeva et al. 2024). Moreover, a wide lensing cluster survey mitigates cosmic variance by sampling random lines of sight, which is a crucial advantage in the search for rare galaxy populations. Importantly, this approach allows us to probe the intrinsically faintest regimes while the detected sources remain relatively bright in the observed frame, significantly accelerating follow-up spectroscopy.

Third, incorporating medium-band filters in JWST legacy NIRCam fields is highly beneficial for ensuring robust Pop III selection. As demonstrated in Figures 1 and 2, F410M is an essential filter for anchoring the underlying continuum and quantifying the Balmer jump as well as the strengths of $H\alpha$ and [O III] lines in conjunction with other broadband filters at $z \simeq 6-7$. While the color–color diagrams presented in Figure 2 rely primarily on broadband filters, apart from F410M, the addition of other medium-band filters further strengthens constraints on $H\alpha$ and/or [O III] line strengths. For example, Appendix B presents color–color diagrams using F480M in place of F444W, demonstrating a similar capability for separating metal-enriched galaxies from Pop III galaxies. A disadvantage with the medium-band filters is their narrower validated redshift range for capturing specific emission lines compared to broadband filters. Nevertheless, the NIRCam colors are sensitive to noise fluctuation, which is particularly significant when exploring faint objects. Therefore, the independent measurements of key emission-line strengths provided by medium-band filters, in addition to the broadband filters, enhance the reliability of identifications in these challenging regimes.

Finally, it is important to emphasize that the nondetection of Pop III galaxies at $z \simeq 6-7$ will also be highly constraining, even as deep NIRCam surveys and follow-up spectroscopy continue to expand. In Figure 11, we find that the volume densities of our Pop III galaxy candidates are consistent with several simulation predictions. This implies that replacing detections with upper limits, which could become significantly more stringent in future surveys, will provide critical constraints on Pop III formation and evolution models. These constraints—where none exist at the moment—will revolutionize our understanding of the nature of the first galaxies.

## 7. Summary

In this paper, we present a novel and efficient NIRCam-based selection method for identifying Pop III galaxies at $z \simeq 6-7$, validated through dedicated completeness and contamination simulations. Leveraging deep NIRCam imaging data from a total area of $\simeq 500 \text{ arcmin}^2$ across several JWST legacy fields, including GLIMPSE, UNCOVER, CEERS,

PRIMER, and JOF, we systematically searched for Pop III galaxy candidates. Below, we summarize the main findings:

1. We developed a fiducial selection method combining color-based and SED-based approaches. The color-based method incorporates three color–color diagrams, optimized to capture the key features of Pop III galaxies at $z = 5.6–6.6$, such as the absence of detectable metal lines like [O III], strong $H\alpha$ emission, and a prominent Balmer jump. The SED-based method relies on the $\chi^2$ difference between best-fit SEDs using Pop III and metal-enriched galaxy templates. Monte Carlo simulations were performed to evaluate completeness and contamination rates. By combining the color and SED-based methods, we confirmed high completeness ($>50\%–80\%$) and zero contamination across most of the parameter space at $z = 5.6–6.6$, above the detection limit of the data (Figures 1–4; Section 2).

2. Under realistic assumptions, we demonstrate that only the deepest JWST imaging surveys are able to breach the luminosity regime where Pop III-dominated galaxies are most likely to occur ($> 30.5$ mag), underscoring the importance of gravitational lensing for meaningful Pop III galaxy searches (Figure 5; Section 2.6).

3. We identified one robust Pop III galaxy candidate, GLIMPSE-16043, at $z = 6.50^{+0.03}_{-0.24}$, and one tentative candidate, JOF-21739, at $z = 6.17^{+0.19}_{-0.06}$. GLIMPSE-16043 satisfies all fiducial selection criteria, whereas JOF-21739 meets the SED-based thresholds but falls near the borders of the color–color diagram criteria, failing to meet all three (Figures 7, 8, 17, and 18; Section 4 and Appendix D).

4. GLIMPSE-16043 is a moderately lensed galaxy ($\mu = 2.9^{+0.1}_{-0.2}$), exhibiting textbook features expected of a Pop III galaxy: strong $H\alpha$ emission with observed rest-frame equivalent width of $2810 \pm 550 \AA$, a pronounced Balmer jump, no detectable metal lines ([O III]/H$\beta < 0.44$), and nascent stellar populations with $t_{age} < 5$ Myr. The inferred stellar mass is $\sim 10^5 M_\odot$, with intrinsic UV magnitudes of $M_{UV} = -15.89^{+0.12}_{-0.14}$. It has an extremely compact morphology consistent with a point source, with an upper limit on its effective radius of $r_e < 40$ pc after the lensing correction (Figures 9 and 19; Section 4.2 and Appendix D).

5. Using the $H\alpha$ and [O III]+H$\beta$ flux excesses observed in F480M and F356W, we constrain the gas-phase metallicity to $12 + \log(O/H) < 6.4$ ($\simeq 0.005 Z_\odot$). In the $M_{UV}-Z_{gas}$ relation, GLIMPSE-16043 significantly deviates from the extrapolated trends of JWST measurements at $z = 4–10$, appearing $>10 \times$ brighter than expected for its metallicity. This deviation may indicate a top-heavy IMF, pristine gas inflows, or extremely metal-poor AGNs (Figure 10, Section 4.2.4).

6. The tentative candidate, JOF-21739 at $z = 6.17^{+0.19}_{-0.06}$ found in the JOF field, displays key SED signatures of Pop III galaxies similar to GLIMPSE-16043: faint ($M_{UV} = -17.62^{+0.17}_{-0.15}$), very young ($t_{age} < 5$ Myr), point-source morphology ($r_e < 80$ pc), and very low metallicity ($12 +$, i.e., $Z_{gas} < 0.003 Z_\odot$) inferred from the F356W and F444W flux excesses. However, the Prospector fitting analysis yields only a small $\Delta\chi^2 = 1.2$ between a $z = 6.2$ Pop III and a $z \sim 1.5$
Paa solution. This reaffirms the tentative nature of JOF-21739, which does not perfectly pass our color–color thresholds (Figures 6, 17, 18, and 19; Appendix D).

7. We derived the first observational constraint on the high-redshift Pop III UVLF and SFRD at $z = 5.7–6.6$. The volume density and SFRD of our best candidate, GLIMPSE-16043, align well with theoretical predictions, providing independent support for its robustness as a Pop III galaxy candidate (Figures 11 and 12; Sections 4.3 and 4.4).

8. While alternative scenarios such as extremely metal-poor galaxies or AGNs cannot be completely ruled out at the moment, even these interpretations make GLIMPSE-16043 and JOF-21739 exceptional laboratories for studying the first stages of galaxy or black hole formation in the early Universe. The low-$z$ scenario of a Paschen-$\alpha$ interloper at $z \approx 1.5$ merits caution due to a nearby neighbor at the same redshift, but this would have to be a hitherto unobserved kind of galaxy that is yet to be discovered. Deep, high-resolution follow-up spectroscopy will be critical for confirming the Pop III nature of these candidates (Figure 14; Section 5).

This work paves a clear path for the discovery of the first Pop III galaxies. Whatever the fate of the present candidates, the methods developed in this study will empower Pop III galaxy searches throughout the JWST era. The volume density ($\approx 10^{-4}$ cMpc$^{-3}$) implied by the finding of our robust candidate is remarkably high, suggesting more than one Pop III source awaiting discovery per NIRCam pointing. The challenge, however, lies in pushing past $\approx 30–31$ mag efficiently, while also remaining sensitive to the Pop III signatures.

Inspired by these results, we project that a snapshot NIRCam survey of the most magnified regions across a large number of lensing clusters, and deploying carefully chosen medium bands, could provide significant observational constraints to help distinguish between different Pop III galaxy simulations, and learn more about the key physical mechanisms of these earliest galaxies and their Pop III star formation.

Exactly 100 yr ago, our cosmic horizon expanded past the edges of the Milky Way for the first time, with Andromeda and Triangulum marking the boundaries of our place in the Universe (e.g., E. P. Hubble 1925). As we reflect on the profound discoveries of the last 100 yr, it is intriguing to consider how those early surveyors of glass plates would view the prospect that we may soon detect the Universe’s very first stars.

## Acknowledgments

We are grateful to the CEERS, PRIMER, JOF, UNCOVER, and GLIMPSE teams for developing their NIRCam surveys, and to the various JWST and HST surveys acknowledged in Section 3 that enabled our search. We thank Kimihiko Nakajima and Kohei Inayoshi for sharing Pop III and/or AGN templates, Steven Finkelstein for comments on the completeness and contamination rate simulation, Aaron Yung for SEDs of simulated galaxies, Joel Leja, Ben Johnson, and Sandro Tacchella for advise on SED fitting, and Takashi Kojima and Hiroto Yanagisawa for discussions.

We made extensive use of the DAWN JWST Archive for various comparisons presented in this paper. Some of the data

products presented herein were retrieved from the Dawn JWST Archive (DJA). DJA is an initiative of the Cosmic Dawn Center (DAWN), which is funded by the Danish National Research Foundation under grant DNRF140. The prism spectra used in this paper were observed as part of the following programs, and we are grateful to these teams for helping build the rich spectroscopic legacy of JWST: 1180, 1181, 1210, 1286, 3215 (A. J. Bunker et al. 2024; F. D'Eugenio et al. 2024); 1211–1215 (M. V. Maseda et al. 2024); 1345 (S. L. Finkelstein et al. 2024); 1433 (T. Y.-Y. Hsiao et al. 2024); 1747 (G. Roberts-Borsani et al. 2025); 2028 (X. Wang et al. 2024); 2073 (PI: J. Hennawi); 2198 (L. Barrufet et al. 2025); 2282 (L. D. Bradley et al. 2023); 2561 (R. Bezanson et al. 2024; S. H. Price et al. 2024); 2565 (T. Nanayakkara et al. 2023); 2750 (P. Arrabal Haro et al. 2023b); 2756 (PI: W. Chen); 2767 (C. C. Williams et al. 2023); 3073 (M. Castellano et al. 2024); 4106 (PI: E. Nelson); 4233 (A. de Graaff et al. 2025); 4446 (B. L. Frye et al. 2024); 4557 (PI: H. Yan); 6541 (PI: E. Egami); 6585 (PI: D. Coulter).

This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. The specific observations can be accessed via doi: 10.17909/xpxt-a441. These observations include data associated with program No. 03293. Support for program No. 03293 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127.

This project has received funding from NASA through the NASA Hubble Fellowship grant HST-HF2-51505.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. This work has received funding from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract No. MB22.00072, as well as from the Swiss National Science Foundation (SNSF) through project grant 200020_207349. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. H.A. and I.C. acknowledge support from CNES, focused on the JWST mission, and the Programme National Cosmology and Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. I.C. acknowledges funding support from the Initiative Physique des Infinis (IPI), a research training program of the Idex SUPER at Sorbonne Université. A.Z. acknowledges support by grant No. 2020750 from the United States–Israel Binational Science Foundation (BSF) and grant No. 2109066 from the United States National Science Foundation (NSF); and by the Israel Science Foundation grant No. 864/23. P.N. acknowledges support from the Gordon and Betty Moore Foundation and the John Templeton Foundation that fund the black hole Initiative (BHI) at Harvard University, where she serves as one of the PIs. B.L. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1—390900948 (the Heidelberg STRUCTURES Excellence Cluster). Y.S. and G.M. have received funding from the European
Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 833925, project STAREX).

Facility: JWST.

Software: SOURCE EXTRACTOR (E. Bertin & S. Arno-uts 1996), ASTROPY (Astropy Collaboration et al. 2013).

## Appendix A
### SED Models

To develop the color-color diagrams for the color-based selection (Section 2.1), we use several Pop III SED models in the literature and generate Pop II SED models using a public SED code. Below, we briefly describe the models we use in our analysis.

#### A.1. Pop III Models

##### A.1.1. E. Zackrisson et al. (2011)

The stellar component is based on single stellar population of Pop III stars from D. Schaefer (2002), D. Schaefer (2003), and A. Raiter et al. (2010) with a total mass of $10^6 M_\odot$, while the nebular emission is computed using the photoionization code Cloudy (G. J. Ferland et al. 1998, 2013, 2017). Yggdrasil supports three types of IMFs for Pop III stars: a very top-heavy IMF using the D. Schaefer (2002) stellar SSP with a power-law IMF with a slope of $\alpha = 2.35$ across the mass range $50-500 M_\odot$; a moderately top-heavy IMF assuming a log-normal distribution with characteristic mass at $10 M_\odot$ and $1 M_\odot$, but with wings extending from $1-500 M_\odot$; and the Kroupa IMF covering a broad stellar mass range of $0.1-100 M_\odot$. These three models are dubbed Pop III.1, Pop III.2, and Pop III. Kroupa, generally reflecting the characteristic mass of stars $M_{\star,IMF}$ to be ~100, 10, and $1 M_\odot$, respectively. The templates span a wide range of stellar ages ($t_{age}$) from 0.01-100 Myr. For this study, we specifically adopt ages of [0.01, 0.1, 1, 5, 10, 30] Myr when available for all three IMF scenarios, with the gas covering fraction of $f_{cov} = 1.0$, ensuring maximal contributions from nebular emission. The Cloudy calculation is performed with the electron density of $n_e = 100 \text{ cm}^{-3}$ in the Yggdrasil model.

##### A.1.2. K. Nakajima & R. Maiolino (2022)

With Cloudy-based calculations, K. Nakajima & R. Maiolino (2022) explored a range of physical conditions for Pop III stellar populations and their surrounding gas, using the SEDs from D. Schaefer (2003) and A. Raiter et al. (2010). These models assume a Salpeter IMF, while varying the mass range: $1-100 M_\odot$, $1-500 M_\odot$, and $50-500 M_\odot$. In this paper, we refer to these three models as Pop III.Sal1-100, Pop III. Sal50-500, and Pop III.Sal1-500, respectively. The wider mass ranges, particularly the top-heavy $50-500 M_\odot$ IMF, represent environments where high-mass stars dominate, as expected in metal-free conditions with limited cooling mechanisms. For our analysis, we use models with a gas electron density of $n_e = 10^3 \text{ cm}^{-3}$ and an ionization parameter of $\log U = -1.5$. These conditions reflect typical values inferred in star-forming regions at high redshift (e.g., Y. Isobe et al. 2023; N. A. Reddy et al. 2023; R. L. Sanders et al. 2023), ensuring consistency with observed properties of high-energy ionizing sources.

#### A.2. Pop II Models

##### A.2.1. BAGPIPES

The mock SEDs are generated with the default setup of BAGPIPES, which assumes a P. Kroupa & C. M. Boily (2002) IMF and the G. Bruzual & S. Charlot (2003) stellar population library. For this study, we explore a parameter space motivated by recent JWST observations of early galaxies, covering gas-phase metallicities of $Z_{gas}/Z_\odot = [0.01, 0.05, 0.10, 0.20]$, ionization parameters $\log U$ of [-3, -2, -1], a fixed electron density of $n_e = 10^3 \text{ cm}^{-3}$, zero dust attenuation, and redshifts $z = 5.0-7.5$ in steps of $\Delta z = 0.1$. These conditions are also optimal to encompass the diversity of metal-enriched young galaxies in the early Universe, which may have NIRCam colors close to the Pop III galaxies. Photometric data for the JWST/NIRCam filters are generated directly from these SEDs, including contributions from both stellar continuum and nebular emission lines.

## Appendix B
### Color-Color Selection with F480M

As discussed in Section 2.1, the selection criteria for Pop III galaxies primarily rely on color-color diagrams constructed using NIRCam broadband filters (F200W, F277W, F356W, and F444W) and the medium-band filter F410M. However, other medium-band filters such as F480M can also play a critical role in refining the selection function by providing additional constraints on the key emission line of $H\alpha$ at similar redshifts, although the available redshift range is even narrowed naturally due to its narrower filter response, compared to F444W.

In Figure 15, we present color-color diagrams incorporating F480M, replacing F444W. We define the color selection enclosed by the following vertices in the color-color space:

1. F277W - F356W versus F410M - F480M (vertices):

$$ \begin{pmatrix} -1.5, 4.5 \\ -1.5, -0.2 \\ -0.1, 1.9 \\ 0.3, 1.9 \\ 0.7, 4.5 \end{pmatrix}, \begin{pmatrix} -0.15, -0.2 \\ -0.1, 1.9 \\ 0.3, 1.9 \\ 0.7, 4.5 \end{pmatrix} \quad (B1) $$

2. F200W - F410M versus F356W - F480M (vertices):

$$ \begin{pmatrix} -2.5, 3.2 \\ -2.5, -0.1 \\ -0.7, 1.5 \\ 0.3, 1.5 \\ 0.5, 3.2 \end{pmatrix}, \begin{pmatrix} -0.8, -0.1 \\ -0.7, 1.5 \\ 0.3, 1.5 \\ 0.5, 3.2 \end{pmatrix} \quad (B2) $$

3. F356W - F410M versus F410M - F480M (vertices):

$$ \begin{pmatrix} 0.5, 3.5 \\ 0.5, 1.5 \\ -2.0, 2.3 \\ -2.0, 3.5 \end{pmatrix}, \quad (B3) $$

The x-axis colors remain sensitive to the Balmer jump and weak [O III] emission, while the y-axis colors are designed to isolate strong $H\alpha$ emission. We find that GLIMPSE-16043 also meets all three color-color diagrams, strengthening its robustness as the Pop III candidate. These diagrams demonstrate that F480M also effectively captures the unique SED features of Pop III galaxies, providing an additional support for distinguishing them from metal-enriched galaxies.
<figure><img src="image_15.png" /><figcaption>Figure 15. Same as Figure 2 (top panels) and Figure 6 (bottom panels), but replaced F444W with F480M, which also efficiently isolates the extreme Hα emission for Pop III galaxies, but a narrower redshift range of z = 6.1–6.6. Among the NIRCam sources at z<sub>phot</sub> = 5.0–7.5 in GLIMPSE (gray dots), no sources meet the three color–color diagrams with F480M, except for GLIMPSE-16043.</figcaption></figure>

## Appendix C
SED-based Selection at z = 3–8

In Section 2.2, we introduced the SED-based method optimized for $z \sim 6–7$, aligning it with the color-based method available at $z = 5.6–6.6$. The strength of the SED-based method, however, is its ability to fully leverage all available photometric data, making it particularly valuable for broader redshift ranges. Notably, the contamination rate in the SED-based selection may be reduced as the number of photometric data points increases. To explore this potential, we extend the SED-based selection to a wider redshift range of $z = 3–8$, assuming a case that all 20 NIRCam broadband and medium-band filters are available, as in surveys like UNCOVER+MegaScience (R. Bezanson et al. 2024; K. A. Suess et al. 2024).

In Figure 16, we present the contamination rate of the SED-based method using mock galaxies from the Santa Cruz Semi-Aanalytic Model catalog (L. Y. A. Yung et al. 2022), in the same manner as Section 2.4, but now extended to the broader redshift range and full filter set. Compared to the filter set in Section 2.2, the contamination rate generally decreases due to the increased photometric coverage. Nevertheless, we observe that most of the parameter space still exhibits contamination rates exceeding 1%, which poses challenges for identifying rare populations like Pop III galaxies. Interestingly, there exists a favorable parameter space at $z \simeq 3–5$ and magnitudes $\lesssim 29.0$, where the contamination rate shows almost entirely

<figure><img src="image_16.png" /><figcaption>Figure 16. Same as the right panel of Figure 3, but for the SED-based method extended to a broad redshift range ($z = 3–8$) using all 20 NIRCam broadband and medium-band filters. Compared to the contamination rate of the SED-based method shown in Figure 3, the inclusion of additional photometric data points reduces the contamination rate across much of the parameter space. However, contamination rates exceeding 1% persist in most regions, presenting challenges for identifying rare populations such as Pop III galaxies.</figcaption></figure>

zero. However, careful simulation studies are necessary to verify whether Pop III galaxies brighter than 29 mag can plausibly exist at $z \simeq 3–5$ (e.g., B. Liu & V. Bromm 2020), making this range a promising but uncertain avenue for future exploration.
## Appendix D
A Tentative Candidate Found in JOF

As introduced in Section 4.1, we find a tentative candidate, JOF-21739, in the JOF field (B. Robertson et al. 2024), which meets the SED criteria, while the NIRCam colors are marginal around the thresholds. Figures 17 and 18 show the NIRCam and HST cutouts and the best-fit SEDs with the Pop III and metal-enriched galaxy templates using *EAZY* for JOF-21739. The faintness in F410M is confirmed also in JOF-21739, supporting the presence of Pop III-like key SED features (Figure 1). JOF-21739 is also detected as a very UV-faint source with $M_{UV} = -17.62_{-0.15}^{+0.17}$ and is well reproduced by a combination of young (0.01–50 Myr) Pop III.1 and Pop III.2 models with no dust attenuation, resulting in $\chi^2 = 8.7$ and $\Delta\chi^2 = 14.5$ compared to the best-fit metal-enriched galaxy SED. Note that, although the $\chi^2$ and subsequent $\Delta\chi^2$ values are generally larger for JOF-21739 compared to GLIMPSE-16043, this difference reflects the availability of more NIRCam filters in the JOF field than in GLIMPSE, used in the SED fitting.

When we run other flexible SED fitting tools of *BAGPIPES*, *BEAGLE*, and *Prospector*, for JOF-21739 in the same

manner as GLIMPSE-16043, we obtain $\Delta\chi^2 = 1.2$ with a low-redshift solution at $z \sim 1.5$ with *Prospector*. In fact, Figure 14 shows that the NIRCam color locations of JOF-21739 overlap with the color distributions of spec-z-confirmed $z \sim 1.5$ galaxies, supporting that the low-z solution serves as a reasonable alternative solution for this candidate. These results reaffirm the less robust nature of JOF-21739 as a Pop III candidate.

Nevertheless, we infer the basic physical properties of JOF-21739 as a Pop III candidate in the same manner as GLIMPSE-16043 as follows. We estimate an H$\alpha$ EW of 3600 ± 430 Å, a mean $t_{age}$ of 0.3 Myr, $\beta = -2.79 \pm 0.05$, and a stellar mass of $\approx 10^{5-6} M_{\odot}$. JOF-21739 also exhibits a compact morphology, and the Sérsic profile fitting with *GALFIT* on the high-resolution F150W image does not converge, reaching the smallest size in the fitting grid. In Figure 19, we show the F150W $1'' \times 1''$ cutout around JOF-21739, together with the PSF and the residual maps, demonstrating that the morphology agrees well with the PSF. We place an upper limit of $\sim 80$ pc for the effective radius, in the same manner as GLIMPSE-16043 (Section 4.2.3). The NIRCam coordinates, photometry, difference in SED codes, and physical properties are summarized in Tables 8, 9, and 10, respectively.

<figure><img src="image_17.png" /><figcaption>**Figure 17.** Same as Figure 7, but for the tentative candidate JOF-21739. The faint F410M is also confirmed in JOF-21739, strongly supporting the presence of a Balmer jump and extremely strong H$\alpha$ (EW $\approx 3600$ Å) that are consistent with Pop III SEDs but are challenging to reproduce along with the other emission-line excesses in metal-enriched galaxy SEDs.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 18.** Same as Figure 8, but for the tentative candidate JOF-21739. The minimal $\Delta\chi^2$ of ~1 obtained in Prospector reaffirms the tentative nature of this candidate.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 19.** Same as Figure 9, but for the tentative candidate of JOF-21739, whose morphology is also consistent with the PSF.</figcaption></figure>

<table><caption>**Table 8**
NIRCam Photometry of JOF-21739</caption><tr><td>R.A. (deg)</td><td>53.0346126</td></tr><tr><td>Decl. (deg)</td><td>-27.893660</td></tr><tr><td>F090W (nJy)</td><td>3.76 ± 1.50</td></tr><tr><td>F115W (nJy)</td><td>10.40 ± 1.22</td></tr><tr><td>F150W (nJy)</td><td>8.22 ± 1.21</td></tr><tr><td>F182M (nJy)</td><td>6.07 ± 0.93</td></tr><tr><td>F200W (nJy)</td><td>6.75 ± 1.11</td></tr><tr><td>F210M (nJy)</td><td>6.21 ± 1.13</td></tr><tr><td>F250M (nJy)</td><td>4.37 ± 1.28</td></tr><tr><td>F277W (nJy)</td><td>6.04 ± 0.81</td></tr><tr><td>F300M (nJy)</td><td>4.86 ± 0.80</td></tr><tr><td>F335M (nJy)</td><td>6.69 ± 0.67</td></tr><tr><td>F356W (nJy)</td><td>6.00 ± 0.79</td></tr><tr><td>F410M (nJy)</td><td>3.54 ± 1.27</td></tr><tr><td>F444W (nJy)</td><td>4.85 ± 0.67</td></tr></table>
<table><caption>**Table 9**
SED Fitting Results for JOF-21739 with Different Codes</caption><thead><tr><th>Fitting Type</th><th>Pop III (EAZY)</th><th>Galaxy (EAZY)</th><th>Galaxy (BEAGLE)</th><th>Galaxy (Prospector)</th><th>Galaxy (BAGPIPES)</th></tr></thead><tbody><tr><td>z<sub>phot</sub></td><td>6.17<sup>+0.06</sup><sub>-0.19</sub></td><td>6.17<sup>+0.19</sup><sub>-0.06</sub></td><td>6.41<sup>+0.14</sup><sub>-0.18</sub></td><td>1.51<sup>+0.03</sup><sub>-0.02</sub></td><td>6.37<sup>+0.16</sup><sub>-0.18</sub></td></tr><tr><td>χ<sup>2</sup></td><td>8.7</td><td>23.2</td><td>13.2</td><td>9.9</td><td>15.8</td></tr></tbody></table>

Note. Although JOF-21739 meets our condition of $\Delta\chi^2 > 9$ with EAZY, its NIRCam colors do not meet all of the color-selection criteria; thereby, we classify it as a tentative Pop III candidate.

<table><caption>**Table 10**
Physical properties of JOF-21739</caption><thead><tr><th>ID</th><th>JOF-21739</th></tr></thead><tbody><tr><td>z<sub>phot</sub></td><td>6.17<sup>+0.19</sup><sub>-0.06</sub></td></tr><tr><td>β</td><td>-2.79 ± 0.05</td></tr><tr><td>M<sub>UV</sub> [mag]</td><td>-17.62<sup>+0.17</sup><sub>-0.15</sub></td></tr><tr><td>M<sub>star</sub> [M<sub>☉</sub>]</td><td>≈10<sup>5</sup>-<sup>6</sup></td></tr><tr><td>r<sub>e</sub> [pc]</td><td>&lt;80</td></tr><tr><td>t<sub>age</sub> [Myr]</td><td>0.3</td></tr><tr><td>H<sub>α</sub> EW [Å]</td><td>3600 ± 430</td></tr><tr><td>O III/Hβ</td><td>&lt;0.32</td></tr><tr><td>12 + log(O/H)</td><td>&lt;6.2</td></tr></tbody></table>

Notes. Physical parameters are derived using the best-fit Pop III templates (see Section 4.2 for details). The lens correction is applied.

ª The $1\sigma$ upper limit is presented, the same as E. Vanzella et al. (2023b).

## Appendix E
Theoretical Frameworks for Pop III UVLF

### E.1. E. Visbal et al. (2020) Model

We utilized a modified version of the semianalytic model from E. Visbal et al. (2020) to estimate the Pop III UVFL. This model uses halo merger trees from cosmological N-body simulations to follow the formation of Pop III and metal-enriched star formation including the three-dimensional variations in Lyman-Werner feedback, reionization, and metal enrichment of the intergalactic medium (IGM). The main change from E. Visbal et al. (2020) is that the minimum halo mass in which Pop III star formation occurs is now calibrated from the simulations of M. Kulkarni et al. (2021). The fitting formulae from M. Kulkarni et al. (2021) include the simultaneous impact of the Lyman-Werner intensity, redshift, and the baryon-dark matter streaming velocity (D. Tseliakhovich & C. Hirata 2010). The model is applied to 10 different 3 Mpc across simulation boxes with dark matter particle resolution of 8000 M☉. This resolution enables us to follow early metal enrichment from Pop III star formation in low-mass mini-halos.

In order to estimate the Pop III UVLF, we determine the number of Pop III galaxies that form in our simulation boxes from $z = 5.9-6.8$ (Pop III star formation is assumed to occur in instantaneous bursts). We convert Pop III stellar mass to UV luminosity with the results from A. Raiter et al. (2010) and D. Schaerer (2002; assuming the IMF in “Model C”). Additionally, we make the simplifying assumption that the Pop III galaxies are bright for a lifetime of 3 Myr and then abruptly become very faint. We note that our 10 simulation boxes are assigned a range of baryon-dark matter streaming velocities (from 0–3 times the typical value of 30 km s⁻¹ at recombination), with the relative abundance of various relative velocities taken into account.

We compare two model parameterizations to the observational constraints in Figure 11. We note that in both models, Lyman-Werner feedback results in all Pop III galaxies found forming in halo masses near the atomic cooling threshold, which are forming stars for the first time. For the first model, we assume the fiducial parameters given in E. Visbal et al. (2020; see Table 1). Here, the star formation efficiency $\epsilon_{*,III}$ (defined as the fraction of gas in a halo that forms Pop III stars during the burst) is 0.001. This leads to Pop III stellar masses of $\sim 10^4$ M☉, which are substantially fainter than GLIMPSE-16043. For the second model, we adopt the same fiducial parameters except for the star formation efficiency, which is increased to 0.01. This leads to Pop III galaxies with stellar mass of $\sim 10^5$ M☉ and a Pop III UVLF that is in remarkable agreement with the constraints from GLIMPSE-16043. We note that from the E. Visbal et al. (2020) model, we expect a second population of Pop III galaxies with luminosities similar to JOF-21739 for an efficiency of 0.01; however, the effective volume of our 10 N-body simulations is not sufficient to probe their abundance. This second population is expected to form in halos more massive than the atomic cooling limit at a virial mass where gas photoheated by reionization can collapse gravitationally.

### E.2. B. Liu & V. Bromm (2020) Model

Starting from the rate density of Pop III star-forming halos at $z \simeq 6.5$, as found in the simulations of B. Liu & V. Bromm (2020), we derive approximate values for the Pop III UV luminosity function for sources with stellar mass $\sim 10^5$ and $\sim 10^4$ M☉, as follows. The emergence rate of Pop III star-forming halos is regulated by the photoheating feedback of reionization and turbulent metal mixing. The latter can be modeled as a diffusion process with coefficient $D_{\text{mix}} = \beta_{\text{mix}}(v_{\text{vir}}R_{\text{vir}}/3)$, where $\beta_{\text{mix}} \lesssim 1$ is an efficiency parameter, and $v_{\text{vir}}$, $R_{\text{vir}}$ are the virial velocity and radius of the host halo, where the Pop III pocket is located. Without enhanced photoheating feedback and metal mixing, corresponding to late reionization (or strong shielding) and $\beta_{\text{mix}} = 0$ (original simulation results), the simulations predict a rate density of $\dot{n} \sim 6 \times 10^{-8} \text{ yr}^{-1} \text{cMpc}^{-3}$. If reionization is nearly complete at $z \simeq 6.5$ (preventing star formation in small halos $\lesssim 4 \times 10^8$ M☉) and metal mixing is highly efficient (with $\beta_{\text{mix}} = 1$), this value can be reduced by up to a factor of $\sim 100$ (see Figure 13 in B. Liu & V. Bromm 2020).

If Pop III clusters are visible for $t \sim 3$ Myr (see Figure 5), the corresponding number density of Pop III-hosting halos is $\dot{n}t \sim 0.2 \text{ cMpc}^{-3}$. The question now is which fraction of these halos host massive Pop III clusters with $\sim 10^5$ M☉ (which is the Pop III stellar mass inferred for GLIMPSE-16043, corresponding to $M_{\text{UV}} \simeq -15.9$, according to the analysis above). It is challenging to answer this question because the formation efficiency and mass distribution of Pop III star clusters are determined by the subparsec-scale properties of Pop III star-forming clouds and the balance between cloud collapse and protostellar feedback (e.g., B. Liu et al. 2024), which are
mostly unresolved in cosmological simulations and treated
with subgrid models. In the B. Liu & V. Bromm (2020)
simulations, such massive clusters do not exist, possibly
because of the limited volume of the computational box and
idealized subgrid models for star formation and stellar
feedback, with an average mass of Pop III stars per halo of
~1000 $M_{\odot}$. For simplicity, we fix the overall formation
efficiency of Pop III stars to that found in the simulations
and group all Pop III stars into $10^5 M_{\odot}$ clusters to estimate the
fraction effectively as ~0.01. The resulting maximum number
density of $10^5 M_{\odot}$ Pop III clusters is $2 \times 10^{-3} \text{cMpc}^{-3}$, which
we take as an approximation for the upper limit of UVLF $\Phi$ at
$M_{\text{UV}} \simeq -15.9$, assuming a magnitude spread of order unity
(see Figure 5). Similarly, we estimate the $\Phi$-value corresp-
onding to a Pop III stellar mass of $\sim 10^4 M_{\odot}$, by assuming
an effective host halo fraction of ~0.1, and a resulting
UV luminosity that is about 2.5 magnitudes fainter
($M_{\text{UV}} \simeq -13.4$). Given these upper limits and considering
the theoretical uncertainty from photoheating feedback and
metal mixing, we represent the UVLF estimates based on
B. Liu & V. Bromm (2020) with the shaded yellow region in
Figure 11. As can be seen, there is reasonable agreement, also
demonstrating the point that the empirical UVLF can serve to
constrain various physical processes governing Pop III star
formation during reionization.

E.3. Sarmento et al. (2022, 2025) Model

We used a customized version of Ramses-RT
(R. Teyssier 2002; J. Rosdahl et al. 2013; J. Rosdahl &
R. Teyssier 2015), a cosmological adaptive mesh refine-
ment simulation to track Pop III star formation in a
24 Mpc h⁻¹ box. We set the maximum refinement level
such that the best average resolution is 91.6 h⁻¹ comoving
pc (cpc) in the densest, refined regions.

Our version of Ramses-RT tracks the unpolluted fraction of
gas in each simulation cell using a subgrid estimate of turbulence.
This allows our team to track Pop III star formation in otherwise
polluted cells (R. Sarmento et al. 2017) by estimating the fraction
of unpolluted gas as a function of time. This model thereby
improves the estimate of Pop III stellar mass at each epoch.

The rest-frame UV (1500 Å) luminosities of our star
particles' (SPs') are based on a set of simple stellar population
SED models parameterized by the SPs' ages, metallicities, and
masses. The SEDs are based on STARBURST 99 (C. Leitherer
et al. 1999) for Population II (Pop II) SPs. Pop III SPs use
A. Raiter et al. (2010) and D. Schaerer (2003) SEDs. Pop II
stars with $Z_{*} > Z_{\text{crit}} = 10^{-5}Z_{\odot}$ are modeled on a E. E. Salpeter
(1955) IMF with masses between 0.8 and 100 $M_{\odot}$. Pop III SPs
with $Z_{*} \le Z_{\text{crit}}$ are based on the A. Raiter et al. (2010) and
D. Schaerer (2003) SEDs for a zero-metallicity population and
utilize a log-normal IMF centered on a characteristic mass of
$60 M_{\odot}$ with $\sigma = 1.0$ and a mass range $1 M_{\odot} \le M_{*} \le 500 M_{\odot}$.
The SEDs model instantaneously bursts across the age range
of SPs in the simulation for both types of stars. The
methodology is outlined in detail in R. Sarmento et al. (2018).

E.4. Venditti et al. (2023, 2024) Model

We adopted the suite of *dustyGadget* (L. Graziani et al. 2020) cosmological simulations introduced in C. Di Cesare et al. (2023), which consists of eight volumes with a comoving side of $50h^{-1}$ cMpc, a total number of $2 \times 672^3$ particles, and

a mass resolution for dark matter/gas particles of
$3.53 \times 10^7 h^{-1} M_{\odot}/5.56 \times 10^6 h^{-1} M_{\odot}$ each, evolved from
$z \simeq 100$ down to $z \simeq 4$. These are the largest simulations
available that include a model for Pop III star formation and
feedback: Pop III particles representing stellar populations
with a mass of $\sim 2 \times 10^6 M_{\odot}$ are formed when the star-forming
condition is met in a gas below the critical metallicity
$Z_{\text{crit}} = 10^{-4} Z_{\odot}$; a Salpeter-like IMF within the mass range
$[100, 500] M_{\odot}$ is assumed, following mass-dependent yields
from stars in the PISN range $[140, 260] M_{\odot}$. The simulations
have demonstrated good agreement with available model
predictions and observations of the cosmic star-formation-
rate/stellar-mass-density evolution and with important scaling
relations (C. Di Cesare et al. 2023). The main findings on the
statistics and physical properties of Pop III-forming environ-
ments have been detailed in A. Venditti et al. (2023), while
their detectability through the He II line at 1640 Å and through
PISNe has also been discussed in A. Venditti et al.
(2024a, 2024b), respectively.

The UVLF estimates in Figure 11 have been obtained by
associating each Pop III particle with the specific luminosity at
1500 Å of the closest Yggdrasil model, i.e., the closest-age SED
from the instantaneous burst, Pop III.1 database, assuming a
Salpeter-like IMF in the range [50, 500] $M_{\odot}$ and zero metallicity,
and including nebular corrections with a covering fraction
$f_{\text{cov}} = 1$. As the Yggdrasil database only presents photoionization
calculations for $10^6 M_{\odot}$ stellar populations, results have been
scaled by the total stellar mass of each particle. We also take into
account an efficiency factor $\eta_{\text{III}} = 0.01-0.1$ with respect to our
Pop III mass resolution element $M_{\text{III, res}} \sim 2 \times 10^6 M_{\odot}$ (such that
$M_{\text{III}} = \eta_{\text{III}} M_{\text{III, res}}$), as described in A. Venditti et al.
(2024a, 2024b). Note that, while these assumptions are not
entirely consistent with the adopted feedback model, they
represent a reasonable compromise between consistency with
the underlying simulation assumptions and a more realistic
emission model. We also caution that the results of photoioniza-
tion simulations are not easily scalable by stellar mass in principle,
as changes in the normalization of the radiation input into the
nebula may result in unpredictable changes in the output spectral
shape. However, we emphasize that the scope of this paper is to
provide a first-order comparison, while a more in-depth study of
the predicted emission of Pop III-hosting *dustyGadget*
galaxies will be provided in a future paper (A Venditti et al.
2025, in preparation), by performing ad hoc *Cloudy* simulations.

All stellar particles from one of the most star-forming cubes in
the simulations (U12) at $z \simeq 6.7$ have been considered here,
without taking into account any halo classification. Particularly,
they also include objects below the nominal stellar-mass-resolution
threshold $M_* \simeq 10^{7.5} M_\odot$ adopted in A. Venditti et al. (2023),
corresponding to galaxies resolved with less than $\sim 20$ stellar
particles. Pop III-hosting halos above this threshold always host
Pop II stellar populations at the same time (with Pop III stars either
found at the periphery of the central galaxy or in metal-free
satellites), while Pop III-only systems are all in the poorly resolved
regime. Tailored zoomed-in simulations will shed further light on
this halo population.

Appendix F
Metal-poor AGN

As discussed in Section 5.2, the compact morphology and
the deviation in the $M_{\text{UV}}-Z_{\text{gas}}$ relation suggest a possible AGN
origin for our candidates. This scenario is particularly
<figure><img src="image_20.png" /><figcaption>**Figure 20.** Best-fit SED with seed BH models presented in K. Inayoshi et al. (2022b) and K. Nakajima & R. Maiolino (2022). The red symbols represent the NIRCam and HST photometry of GLIMPSE-16043 in the same manner as Figure 7. The large deviations are observed in the rest-frame UV continuum shape, resulting in the significant Δχ² values and suggesting that GLIMPSE-16043 is less likely explained by the typical metal-poor AGN models in literature. However, the deviations might be mitigated by incorporating alternative surrounding gas conditions where nebular continuum emission dominates the SED, and we cannot rule out the possibility of the metal-poor AGN scenario.</figcaption></figure>

intriguing, as our candidates may represent seed or intermediate-mass BHs with $M_{BH} \simeq 10^{4-5} M_{\odot}$, given their ultrafaint luminosities. In Figure 20, we present the best-fit SEDs for GLIMPSE-16043 obtained using EAZY and the metal-poor AGN templates introduced by K. Inayoshi et al. (2022b) and K. Nakajima & R. Maiolino (2022). We find that the best-fit AGN SEDs yield significantly higher χ² values (40–100) and subsequent Δχ² values compared to the best-fit SED with Pop III templates (χ² = 2.7; see Figure 7). The primary discrepancy lies in the continuum shape: the AGN templates exhibit a red continuum, which contrasts with the blue continuum observed in our candidates. This difference might be mitigated by assuming alternative surrounding gas conditions where nebular continuum emission dominates the SED. Therefore, although the current results disfavor the AGN scenario, we cannot entirely rule out the possibility of extremely metal-poor AGN origins for our candidates. Further studies that explore AGN models incorporating nebular-dominated emission will be essential for a comprehensive assessment of this scenario.

## ORCID iDs

Seiji Fujimoto https://orcid.org/0000-0001-7201-5066

Rohan P. Naidu https://orcid.org/0000-0003-3997-5705

John Chisholm https://orcid.org/0000-0002-0302-2577

Hakim Atek https://orcid.org/0000-0002-7570-0824

Ryan Endsley https://orcid.org/0000-0003-4564-2771

Vasily Kokorev https://orcid.org/0000-0002-5588-9156

Lukas J. Furtak https://orcid.org/0000-0001-6278-032X

Richard Pan https://orcid.org/0000-0002-9651-5716

Boyuan Liu https://orcid.org/0000-0002-4966-7450

Volker Bromm https://orcid.org/0000-0003-0212-2979

Alessandra Venditti https://orcid.org/0000-0003-2237-0777

Eli Visbal https://orcid.org/0000-0002-8365-0337

Richard Sarmento https://orcid.org/0000-0002-8013-5970

Andrea Weibel https://orcid.org/0000-0001-8928-4465

Pascal A. Oesch https://orcid.org/0000-0001-5851-6649

Gabriel Brammer https://orcid.org/0000-0003-2680-005X

Daniel Schaerer https://orcid.org/0000-0001-7144-7182

Angela Adamo https://orcid.org/0000-0002-8192-8091

Danielle A. Berg https://orcid.org/0000-0002-4153-053X

Rachel Bezanson https://orcid.org/0000-0001-5063-8254

Rychard Bouwens https://orcid.org/0000-0002-4989-2471

Iryna Chemerynska https://orcid.org/0009-0009-9795-6167

Adélaïde Claeysens https://orcid.org/0000-0001-7940-1816

Miroslava Dessauges-Zavadsky https://orcid.org/0000-0003-0348-2917

Anna Frebel https://orcid.org/0000-0002-2139-7145

Damien Korber https://orcid.org/0000-0002-3897-6856

Ivo Labbe https://orcid.org/0000-0002-2057-5376

Rui Marques-Chaves https://orcid.org/0000-0001-8442-1846

Jorryt Matthee https://orcid.org/0000-0003-2871-127X

Kristen B. W. McQuinn https://orcid.org/0000-0001-5538-2614

Julian B. Muñoz https://orcid.org/0000-0002-8984-0465

Priyamvada Natarajan https://orcid.org/0000-0002-5554-8896

Alberto Saldana-Lopez https://orcid.org/0000-0001-8419-3062

Katherine A. Suess https://orcid.org/0000-0002-1714-1905

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

Adi Zitrin https://orcid.org/0000-0002-0350-4488

## References

Adamo, A., Bradley, L. D., Vanzella, E., et al. 2024, Natur, 632, 513

Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023a, Natur, 622, 707

Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023b, ApJL, 951, L22

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJL, 946, L12

Barrufet, L., Oesch, P., Marques-Chaves, R., et al. 2025, MNRAS, 537, 3453

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

Bezanson, R., Labbe, I., Whitaker, K. E., et al. 2024, ApJ, 974, 92

Bouwens, R. J., Smit, R., Schouws, S., et al. 2022, ApJ, 931, 160

Bradley, L. D., Coe, D., Brammer, G., et al. 2023, ApJ, 955, 13

Brammer, G. 2022, gbrammer/msaexp: Full Working Version with 2D Drizzling and Extraction, v0.3, doi:10.5281/zenodo.7299501

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Bromm, V. 2013, RPPh, 76, 112901

Bromm, V., Kudritzki, R. P., & Loeb, A. 2001, ApJ, 552, 464

Bromm, V., & Loeb, A. 2003, Natur, 425, 812

Bromm, V., & Yoshida, N. 2011, ARA&A, 49, 373

Brown, T. M., Tumlinson, J., Geha, M., et al. 2014, ApJ, 796, 91

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2024, A&A, 690, A288

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

Carnall, A. C., Begley, R., McLeod, D. J., et al. 2023a, MNRAS, 518, L45

Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2023b, Natur, 619, 716

Carniani, S., D'Eugenio, F., Ji, X., et al. 2025, A&A, 696, A87

Carniani, S., Hainline, K., D'Eugenio, F., et al. 2024, Natur, 633, 318

Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31

Castellano, M., Napolitano, L., Fontana, A., et al. 2024, ApJ, 972, 143

Chemerynska, I., Atek, H., Furtak, L. J., et al. 2024, MNRAS, 531, 2615

Chevallard, J., & Charlot, S. 2016, MNRAS, 462, 1415

Chon, S., & Omukai, K. 2025, MNRAS, 539, 2561

Chon, S., Omukai, K., & Schneider, R. 2021, MNRAS, 508, 4175
Claeyssens, A., Adamo, A., Messa, M., et al. 2025, MNRAS, 537, 2535

Davé, R. 2008, MNRAS, 385, 147

de Graaff, A., Setton, D. J., Brammer, G., et al. 2025, NatAs, 9, 280

Dessauges-Zavadsky, M., Marques-Chaves, R., Schaerer, D., et al. 2025, A&A, 693, A17

D'Eugenio, F., Pérez-González, P. G., Maiolino, R., et al. 2024, NatAs,

Di Cesare, C., Graziani, L., Schneider, R., et al. 2023, MNRAS, 519, 4632

Donnan, C. T., Dunlop, J. S., McLure, R. J., McLeod, D. J., & Cullen, F. 2025, MNRAS, 539, 2409

Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011

Donnan, C. T., McLure, R. J., Dunlop, J. S., et al. 2024, MNRAS, 533, 3222

Downes, A. J. B., Peacock, J. A., Savage, A., & Carrie, D. R. 1986, MNRAS, 218, 31

Eisenstein, D. J., Johnson, B. D., Robertson, B., et al. 2023a, arXiv:2310.12340

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

Endsley, R., Chisholm, J., Stark, D. P., Topping, M. W., & Whitler, L. 2025, ApJ, 987, 189

Endsley, R., Stark, D. P., Whitler, L., et al. 2023, MNRAS, 524, 2312

Fardal, M. A., Katz, N., Weinberg, D. H., & Davé, R. 2007, MNRAS, 379, 985

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, RMXA, 53, 385

Fland, G. J., Korista, K. T., Verner, D. A., et al. 1998, PASP, 110, 761

Fland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, RMXA, 49, 137

Finkelstein, S. L., Leung, G. C. K., Bagley, M. B., et al. 2024, ApJL, 969, L2

Finkelstein, S. L., Ryan, R. E. J., Papovich, C., et al. 2015, ApJ, 810, 71

Frebel, A., & Norris, J. E. 2015, ARA&A, 53, 631

Frye, B. L., Pascale, M., Pierel, J., et al. 2024, ApJ, 961, 171

Fu, S. W., Weisz, D. R., Starkenburg, E., et al. 2023, ApJ, 958, 167

Fujimoto, S., Ouchi, M., Kohno, K., et al. 2024a, arXiv:2402.18543

Fujimoto, S., Kohno, K., Ouchi, M., et al. 2024b, ApJS, 275, 36

Fujimoto, S., Wang, B., Weaver, J., et al. 2024c, ApJ, 977, 250

Furtak, L. J., Zitrin, A., Weaver, J., et al. 2023, MNRAS, 523, 4568

Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001

Gehrels, N. 1986, ApJ, 303, 336

Graziani, L., Schneider, R., Ginolfi, M., et al. 2020, MNRAS, 494, 1071

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

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

Harikane, Y., Ouchi, M., Oguri, M., et al. 2023a, ApJS, 265, 5

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

Hegde, S., & Furlanetto, S. R. 2023, MNRAS, 525, 428

Heintz, K. E., Brammer, G. B., Watson, D., et al. 2025, A&A, 693, A60

Helton, J. M., Rieke, G. H., Alberts, S., et al. 2025, NatAs, 9, 729

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

Hsiao, T. Y.-Y., Álvarez-Márquez, J., Coe, D., et al. 2024, ApJ, 973, 81

Hubble, E. P. 1925, PA, 33, 252

Inayoshi, K., Nakatani, R., Toyouchi, D., et al. 2022a, ApJ, 927, 237

Inayoshi, K., Onoue, M., Sugahara, Y., Inoue, A. K., & Ho, L. C. 2022b, ApJL, 931, L25

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

Inoue, A. K. 2011, MNRAS, 415, 2920

Isobe, Y., Ouchi, M., Tominaga, N., et al. 2023, ApJ, 959, 100

Izotov, Y. I., Thuan, T. X., & Lipovetsky, V. A. 1997, ApJS, 108, 1

Jaacks, J., Finkelstein, S. L., & Bromm, V. 2019, MNRAS, 488, 2202

Jeon, M., Besla, G., & Bromm, V. 2017, ApJ, 848, 85

Jeon, M., Bromm, V., Besla, G., Yoon, J., & Choi, Y. 2021, MNRAS, 502, 1

Jeon, M., Pawlik, A. H., Bromm, V., & Milosavljević, M. 2014, MNRAS, 444, 3288

Jimenez, R., & Haiman, Z. 2006, Natur, 440, 501

Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22

Johnson, J. L., Greif, T. H., & Bromm, V. 2008, MNRAS, 388, 26

Karlsson, T., Bromm, V., & Bland-Hawthorn, J. 2013, RvMP, 85, 809

Karlsson, T., Johnson, J. L., & Bromm, V. 2008, ApJ, 679, 6

Kashino, D., Lilly, S. J., Matthee, J., et al. 2023, ApJ, 950, 66

Katz, H., Kimm, T., Ellis, R. S., Devriendt, J., & Slyz, A. 2023, MNRAS, 524, 351

Klessen, R. S., & Glover, S. C. O. 2023, ARA&A, 61, 65

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

Kojima, T., Ouchi, M., Rauch, M., et al. 2020, ApJ, 898, 142

Kokorev, V., Atek, H., Chisholm, J., et al. 2025, ApJL, 983, L22

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

Kroupa, P., & Boily, C. M. 2002, MNRAS, 336, 1188

Kulkarni, M., Visbal, E., & Bryan, G. L. 2021, ApJ, 917, 40

Labbé, I., Oesch, P. A., Bouwens, R. J., et al. 2013, ApJL, 777, L19

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Liu, B., & Bromm, V. 2020, MNRAS, 497, 2839

Liu, B., Gurian, J., Inayoshi, K., et al. 2024, MNRAS, 534, 290

Liu, B., Sibony, Y., Meynet, G., & Bromm, V. 2025, ApJL, 980, L30

Llerena, M., Amorin, R., Pentericci, L., et al. 2024, A&A, 691, A59

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

Lusso, E., Worseck, G., Hennawi, J. F., et al. 2015, MNRAS, 449, 4204

Lutz, D., Berta, S., Contursi, A., et al. 2016, A&A, 591, A136

Magg, M., Schauer, A. T. P., Klessen, R. S., et al. 2022, ApJ, 929, 119

Maio, U., Ciardi, B., Dolag, K., Tornatore, L., & Khochfar, S. 2010, MNRAS, 407, 1003

Maio, U., Petkova, M., De Lucia, G., & Borgani, S. 2016, MNRAS, 460, 3733

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

Maiolino, R., Uebler, H., Perna, M., et al. 2024b, A&A, 687, 67

Maseda, M. V., de Graaff, A., Franx, M., et al. 2024, A&A, 689, A73

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

Messa, M., Dessauges-Zavadsky, M., Richard, J., et al. 2022, MNRAS, 516, 2420

Morishita, T., Stiavelli, M., Chary, R.-R., et al. 2024, ApJ, 963, 9

Moriya, T. J., Inserra, C., Tanaka, M., et al. 2022, A&A, 666, A157

Mowla, L., Iyer, K., Asada, Y., et al. 2024, Natur, 636, 332

Muñoz, J. B., Qin, Y., Mesinger, A., et al. 2008, MNRAS, 511, 3657

Nagao, T., Sasaki, S. S., Maiolino, R., et al. 2008, ApJ, 680, 100

Naidu, R. P., Matthee, J., Kramarenko, I., et al. 2024, arXiv:2410.01874

Naidu, R. P., Oesch, P. A., Setton, D. J., et al. 2022a, arXiv:2208.02794

Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022b, ApJL, 940, L14

Nakajima, K., & Maiolino, R. 2022, MNRAS, 513, 5134

Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, ApJS, 269, 33

Nakajima, K., Ouchi, M., Xu, Y., et al. 2022, ApJS, 262, 3

Nanayakkara, T., Glazebrook, K., Jacobs, C., et al. 2023, ApJL, 947, L26

Napolitano, L., Castellano, M., Pentericci, L., et al. 2025, A&A, 693, A50

Natarajan, P. 2021, MNRAS, 501, 1413

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

Nishigaki, M., Ouchi, M., Nakajima, K., et al. 2023, ApJ, 952, 11

Nozawa, T., Kozasa, T., Tominaga, N., et al. 2010, ApJ, 713, 356

Oesch, P. A., Brammer, G., Naidu, R. P., et al. 2023, MNRAS, 525, 2864

Oh, S. P., & Haiman, Z. 2002, ApJ, 569, 558

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Ono, Y., Harikane, Y., Ouchi, M., et al. 2023, ApJ, 951, 72

Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Melville, NY: University Science Books)

Page, T. 1966, Sci, 151, 1411

Pallottini, A., Ferrara, A., Gallerani, S., Salvadori, S., & D'Odorico, V. 2014a, MNRAS, 440, 2498

Pallottini, A., Gallerani, S., & Ferrara, A. 2014b, MNRAS, 444, L105

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097

Prgomet, M., Rey, M. P., Andersson, E. P., et al. 2022, MNRAS, 513, 2326

Price, S. H., Bezanson, R., Labbe, I., et al. 2024, ApJ, 982, 51

Raiter, A., Schaerer, D., & Fosbury, R. A. E. 2010, A&A, 523, A64

Reddy, N. A., Topping, M. W., Sanders, R. L., Shapley, A. E., & Brammer, G. 2023, ApJ, 952, 167

Rigby, J. R., Vieira, J. D., Phadke, K. A., et al. 2025, ApJ, 978, 108

Roberts-Borsani, G., Bagley, M., Rojas-Ruiz, S., et al. 2025, ApJ, 983, 18

Roberts-Borsani, G. W., Bouwens, R. J., Oesch, P. A., et al. 2016, ApJ, 823, 143

Robertson, B., Johnson, B. D., Tacchella, S., et al. 2024, ApJ, 970, 31

Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188

Rosdahl, J., & Teyssier, R. 2015, MNRAS, 449, 4380

Salpeter, E. E. 1955, ApJ, 121, 161

Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2023, ApJ, 955, 54

Sarmento, R., & Scannapieco, E. 2022, ApJ, 935, 174

Sarmento, R., Scannapieco, E., & Cohen, S. 2018, ApJ, 854, 75

Sarmento, R., Scannapieco, E., & Pan, L. 2017, ApJ, 834, 23

Savino, A., Weisz, D. R., Skillman, E. D., et al. 2023, ApJ, 956, 86

Scannapieco, C., Tissera, P. B., White, S. D. M., & Springel, V. 2005, MNRAS, 364, 552

Schaerer, D. 2002, A&A, 382, 28

Schaerer, D. 2003, A&A, 397, 527

Schaerer, D., Guibert, J., Marques-Chaves, R., & Martins, F. 2025, A&A, 639, A271

Schauer, A. T. P., Drory, N., & Bromm, V. 2020, ApJ, 904, 145

Schneider, R., Omukai, K., Bianchi, S., & Valiante, R. 2012, MNRAS, 419, 1566

Schneider, R., Omukai, K., Inoue, A. K., & Ferrara, A. 2006, MNRAS, 369, 1437

Schouws, S., Bouwens, R. J., Ormerod, K., et al. 2025, ApJ, 988, 19
Simmonds, C., Tacchella, S., Hainline, K., et al. 2024, MNRAS, 527, 6139

Simon, J. D. 2019, ARA&A, 57, 375

Skinner, D., & Wise, J. H. 2020, MNRAS, 492, 4386

Smit, R., Bouwens, R. J., Franx, M., et al. 2015, ApJ, 801, 122

Sobral, D., Matthee, J., Darvish, B., et al. 2015, ApJ, 808, 139

Stefanon, M., Bouwens, R. J., Illingworth, G. D., et al. 2022, ApJ, 935, 94

Suess, K. A., Weaver, J. R., Price, S. H., et al. 2024, ApJ, 976, 101

Tang, M., Stark, D. P., Chen, Z., et al. 2023, MNRAS, 526, 1657

Teyssier, R. 2002, A&A, 385, 337

Tornatore, L., Ferrara, A., & Schneider, R. 2007, MNRAS, 382, 945

Trussler, J. A. A., Conselice, C. J., Adams, N. J., et al. 2023, MNRAS, 525, 5328

Tseliakhovich, D., & Hirata, C. 2010, PhRvD, 82, 083520

Tumlinson, J. 2006, ApJ, 641, 1

Tumlinson, J., Giroux, M. L., & Shull, J. M. 2001, ApJL, 550, L1

Valentino, F., Brammer, G., Gould, K. M. L., et al. 2023, ApJ, 947, 20

Vanzella, E., Loiacono, F., Bergamini, P., et al. 2023a, A&A, 678, A173

Vanzella, E., Claeyssens, A., Welch, B., et al. 2023b, ApJ, 945, 53

Venditti, A., Bromm, V., Finkelstein, S. L., Graziani, L., & Schneider, R. 2024a, MNRAS, 527, 5102

Venditti, A., Bromm, V., Finkelstein, S. L., et al. 2024b, ApJL, 973, L12

Venditti, A., Graziani, L., Schneider, R., et al. 2023, MNRAS, 522, 3809

Vikaeus, A., Zackrisson, E., Schaerer, D., et al. 2022, MNRAS, 512, 3030

Visbal, E., Bryan, G. L., & Haiman, Z. 2020, ApJ, 897, 95

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

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

Vujeva, L., Steinhardt, C. L., Jespersen, C. K., et al. 2024, ApJ, 974, 23

Wang, X., Cheng, C., Ge, J., et al. 2024, ApJL, 967, L42

Weaver, J. R., Cutler, S. E., Pan, R., et al. 2024, ApJS, 270, 7

Weibel, A., Oesch, P. A., Barrufet, L., et al. 2024, MNRAS, 533, 1808

Weisz, D. R., Dolphin, A. E., Skillman, E. D., et al. 2014, ApJ, 789, 147

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

Xiao, M., Oesch, P. A., Elbaz, D., et al. 2024, Natur, 635, 311

Xu, H., Norman, M. L., O'Shea, B. W., & Wise, J. H. 2016, ApJ, 823, 140

Yajima, H., & Khochfar, S. 2017, MNRAS, 467, L51

Yung, L. Y. A., Somerville, R. S., Ferguson, H. C., et al. 2022, MNRAS, 515, 5416

Zackrisson, E., Rydberg, C.-E., Schaerer, D., Östlin, G., & Tuli, M. 2011, ApJ, 740, 13