# A Glimpse of the New Redshift Frontier through AS1063

Vasily Kokorev, Hakim Atek, John Chisholm, Ryan Endsley, Iryna

Chemerynska, Julian B. Muñoz, Lukas J. Furtak, Richard Pan, Danielle Berg,

Seiji Fujimoto, et al.

► To cite this version:

Vasily Kokorev, Hakim Atek, John Chisholm, Ryan Endsley, Iryna Chemerynska, et al. A Glimpse of the New Redshift Frontier through AS1063. The Astrophysical Journal Letters, 2025, 983, 10.3847/2041-8213/adc458. insu-05101668

HAL Id: insu-05101668

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

Submitted on 8 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# A Glimpse of the New Redshift Frontier through AS1063

Vasily Kokorev<fnref n="1" />, Hakim Atek<fnref n="2" />, John Chisholm<fnref n="1" />, Ryan Endsley<fnref n="1" />, Iryna Chemerynska<fnref n="2" />, Julian B. Muñoz<fnref n="1" />, Lukas J. Furtak<fnref n="3" />, Richard Pan<fnref n="4" />, Danielle Berg<fnref n="1" />, Seiji Fujimoto<fnref n="1" />,<fnref n="5" />, Pascal A. Oesch<fnref n="5" />,<fnref n="6" />, Andrea Weibel<fnref n="6" />, Angela Adamo<fnref n="7" />, Jeremy Blaizot<fnref n="8" />, Rychard Bouwens<fnref n="9" />, Miroslava Dessauges-Zavadsky<fnref n="6" />, Gourav Khullar<fnref n="10" />, Damien Korber<fnref n="6" />, Ilias Goovaerts<fnref n="11" />, Michelle Jecmen<fnref n="1" />, Ivo Labbé<fnref n="12" />, Floriane Leclercq<fnref n="13" />, Rui Marques-Chaves<fnref n="6" />, Charlotte Mason<fnref n="5" />, Kristen B. W. McQuinn<fnref n="11" />,<fnref n="14" />, Rohan Naidu<fnref n="15" />, Priyamvada Natarajan<fnref n="16" />,<fnref n="17" />, Erica Nelson<fnref n="18" />, Joki Rosdahl<fnref n="8" />, Alberto Saldana-Lopez<fnref n="7" />, Daniel Schaerer<fnref n="6" />, Maxime Trebitsch<fnref n="19" />, Marta Volonteri<fnref n="2" />, and Adi Zitrin<fnref n="3" />

## Abstract

We report the discovery of two galaxy candidates at redshifts between $15.7 < z < 16.4$ in James Webb Space Telescope (JWST) observations from the GLIMPSE survey. These robust sources were identified using a combination of Lyman break selection and photometric redshift estimates. The ultradeep NIRCam imaging from GLIMPSE, combined with the strong gravitational lensing of the AS1063 galaxy cluster, allows us to probe an intrinsically fainter population (down to $M_{UV} = -17.0$ mag) than previously achievable. These galaxies have absolute magnitudes ranging from $M_{UV} = -17.0$ to $-17.2$ mag, with blue ($\beta \simeq -2.87$) ultraviolet (UV) continuum slopes, consistent with young, dust-free stellar populations. The number density of these objects, $\log_{10}(\phi/\text[Mpc}^{-3}\text{mag}^{-1}) = -3.47^{+0.13}_{-0.10}$ at $M_{UV} = -17$, is in clear tension with pre-JWST theoretical predictions, extending the overabundance of galaxies from $z \sim 10$ to $z \sim 17$. These results, together with the scarcity of brighter galaxies in other public surveys, suggest a steep decline in the bright end of the UV luminosity function at $z \sim 16$, implying efficient star formation and possibly a close connection to the halo mass function at these redshifts. Testing a variety of star formation histories suggests that these sources are plausible progenitors of the unusually UV-bright galaxies that JWST now routinely uncovers at $z = 10-14$. Overall, our results indicate that the luminosity distribution of the earliest star-forming galaxies could be shifting toward fainter luminosities, implying that future surveys of cosmic dawn will need to explore this faint luminosity regime.

Unified Astronomy Thesaurus concepts: High-redshift galaxies (734); Early universe (435)

## 1. Introduction

According to the standard paradigm of structure formation, the same primordial fluctuations that gave rise to hot and cold spots in the cosmic microwave background will eventually grow, collapse, and form the first galaxies during cosmic dawn, ushering in the epoch of first light (e.g., A. Loeb & S. R. Furlanetto 2013). These first galaxies have remained outside of our observational reach for decades. That is because they are faint and highly redshifted. Even the deepest Hubble Space Telescope (HST) surveys have fallen short of observing first light (P. A. Oesch et al. 2016). The

ultraviolet (UV) light from the first galaxies drops precipitously owing to absorption by foreground neutral gas (the Ly$\alpha$ break; see, e.g., C. C. Steidel et al. 1996), which makes galaxies above $z \gtrsim 10$ invisible in the HST IR filters.

The James Webb Space Telescope (JWST) was designed to observe the first galaxies. With its enormous collecting area and unprecedented near-infrared imaging and spectroscopic capabilities, JWST/NIRCam is sensitive to faint light up to $5 \ \mu\text{m}$ (M. J. Rieke et al. 2023). These are exactly the wavelengths required to discover galaxies forming stars during the first few hundred million years of cosmic history, as their Ly$\alpha$ break falls at $2.0$ and $2.5 \ \mu\text{m}$ at redshifts 16 and 20, respectively. This means that broad imaging bands, such as F200W and F277W, can be paired to search for continuum “dropouts” up to $z \sim 20$, a mere $180$ million years after the big bang.

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA; vkokorev@utexas.edu</footnote>

<footnote n="2">²Institut d'Astrophysique de Paris, CNRS, Sorbonne Université, 98bis Boulevard Arago, 75014, Paris, France</footnote>

<footnote n="3">³Department of Physics, Ben-Gurion University of the Negev, P.O. Box 653, Be'er-Sheva 84105, Israel</footnote>

<footnote n="4">⁴Department of Physics & Astronomy, Tufts University, Medford, MA 02155, USA</footnote>

<footnote n="5">⁵Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen, Jagtvej 128, København N, DK-2200, Denmark</footnote>

<footnote n="6">⁶Département d'Astronomie, Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland</footnote>

<footnote n="7">⁷Department of Astronomy, The Oskar Klein Centre, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden</footnote>

<footnote n="8">⁸Université Claude Bernard Lyon 1, CRAL UMR5574, ENS de Lyon, CNRS, Villeurbanne, F-69622, France</footnote>

<footnote n="9">⁹Leiden Observatory, Leiden University, NL-2300 RA Leiden, Netherlands</footnote>

<footnote n="10">¹⁰Department of Physics & Astronomy and PITT PACC, University of Pittsburgh, Pittsburgh, PA 15260, USA</footnote>

<footnote n="11">¹¹Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA</footnote>

<footnote n="12">¹²Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Melbourne, VIC 3122, Australia</footnote>

<footnote n="13">¹³Univ Lyon, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France</footnote>

<footnote n="14">¹⁴Department of Physics & Astronomy, Rutgers, The State University of New Jersey, Piscataway, NJ 08854, USA</footnote>

<footnote n="15">¹⁵MIT Kavli Institute for Astrophysics and Space Research, 70 Vassar Street, Cambridge, MA 02139, USA</footnote>

<footnote n="16">¹⁶Department of Astronomy, Yale University, 219 Prospect Street, New Haven, CT 06511, USA</footnote>

<footnote n="17">¹⁷Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="18">¹⁸Department for Astrophysical and Planetary Science, University of Colorado, Boulder, CO 80309, USA</footnote>

<footnote n="19">¹⁹LERMA, Sorbonne Université, Observatoire de Paris, PSL Research University, CNRS, 75014 Paris, France</footnote>

<footnote>Received 2024 November 27; revised 2025 March 5; accepted 2025 March 21; published 2025 April 8</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** GLIMPSE field overview. Positions of the high-z candidates are overlaid on top of the red, green, and blue image constructed from all available JWST broad bands. All of the high-z candidates are located in the full-coverage (six primary dithers) area of Module A (left). We additionally overlay a z = 16 critical curve in red.</figcaption></figure>

JWST has successfully used this technique to discover dozens of galaxies at $10 < z < 15$. These observations have revealed a surprising overabundance of UV-bright galaxies at $z > 10$, challenging pre-JWST predictions (to name just a few; N. J. Adams et al. 2023; H. Atek et al. 2023; D. Austin et al. 2023; E. Curtis-Lake et al. 2023; C. T. Donnan et al. 2023; S. L. Finkelstein et al. 2024; P. G. Pérez-González et al. 2023; S. Carniani et al. 2024b; M. Castellano et al. 2024; I. Chemerynska et al. 2024; D. J. McLeod et al. 2024; B. Robertson et al. 2024). Attempts to explain this overabundance of bright galaxies at extreme redshifts invoke extremely efficient star formation in the early Universe, bursty star formation histories (SFHs), a lack of dust attenuation, a top-heavy initial mass function (IMF), more efficient formation of dark matter (DM) halos, or modifications to the ΛCDM paradigm (e.g., F. Pacucci et al. 2022; M. Boylan-Kolchin 2023; A. Dekel et al. 2023; A. Ferrara et al. 2023; S. L. Finkelstein et al. 2024; Y. Harikane et al. 2023; F. Sun et al. 2024; Z. Li et al. 2024).

The ancestors of these UV-bright galaxies at even higher redshift ($z > 15$) have remained elusive. Many observational campaigns have searched for detections of first-light galaxies, but few candidates pass rigorous scrutiny. Extremely dusty star-forming galaxies (R. P. Naidu et al. 2022; P. Arrabal Haro et al. 2023; J. A. Zavala et al. 2023) and spurious detections can masquerade as faint F200W dropout galaxies. A key issue with detecting first-light galaxies is that they are expected to be intrinsically faint, as they have had little time to assemble

significant stellar mass. The discovery and characterization of these elusive primordial galaxies provide empirical constraints on the astrophysics shaping their formation.

Here we report on observations from the JWST GLIMPSE survey (H. Atek et al. 2025, in preparation). GLIMPSE was specifically aimed to detect galaxies during the epoch of first light, by pairing ultradeep NIRCam imaging across seven wide and two medium bands with the gravitational lensing of the foreground galaxy cluster AS1063. With even modest magnification factors of 2, the 30.8 mag GLIMPSE observations can reach absolute magnitudes ($M_{UV}$) of $\sim -17$ at $z \sim 16$, probing the faint galaxies that likely existed during the epoch of first light. We use these observations to identify two robust $z \gtrsim 16$ candidates, observed within a single NIRCam pointing (Figure 1), and assess the validity of these candidates via model fitting, morphology, and properties of their stellar populations.

Throughout this work we assume a flat $\Lambda$CDM cosmology with $\Omega_{m,0} = 0.3$, $\Omega_{\Lambda,0} = 0.7$, and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$ and a G. Chabrier (2003) IMF between 0.1 and 100 $M_{\odot}$. All magnitudes are expressed in the AB system (J. B. Oke 1974).

## 2. Observations and Data

Detailed descriptions of the observations, data reduction, cluster light removal, and source extraction will be presented in the GLIMPSE survey (H. Atek et al. 2025, in preparation) and are briefly summarized below.
## 2.1. GLIMPSE Survey

This work uses the ultradeep imaging from the public GLIMPSE survey (PID: 3293; PIs: H. Atek & J. Chisholm). GLIMPSE targets AS1063, one of the highest-magnification regions in the Hubble Frontier Fields (HFF; J. M. Lotz et al. 2017), with seven broadband filters (F090W, F115W, F150W, F200W, F277W, F356W, F444W) and two medium-band filters (F410M, F480M). The total duration of the survey is ~155 hr of science time, reaching unprecedented observed depths, down to 30.6 mag uniformly across all wide bands. Specifically, the bands used in the $z > 15$ dropout selection, F200W and F277W, have integration times of 19 and 23 hr, respectively. Furthermore, the lensing magnification of AS1063 also means that we can probe intrinsically faint sources that would otherwise be invisible in the deepest JWST blank-field surveys.

GLIMPSE observations use a MEDIUM8 readout pattern for all exposures and adopt a six-position primary dither pattern to cover the short-wavelength (SW) intramodule gaps while maximizing the full-depth area. In addition, we use a subpixel dither with four positions to best sample the point-spread function (PSF). This was done to optimize the signal-to-noise ratio (S/N), with dithers large enough to mitigate the fixed pattern noise, imperative when searching for high-z targets.

## 2.2. Data Reduction

NIRCam imaging data for all seven broadband and two medium-band filters in AS1063 are reduced following the procedure in R. Endsley et al. (2024) using the *jwst_1293_pmap* context map. We implement crucial enhancements over the standard STScI pipeline, including corrections for cosmic rays, stray light, 1/f noise, and detector artifacts (L. D. Bradley et al. 2023; J. Rigby et al. 2023). Furthermore, we construct our own set of image flats based on all public NIRCam imaging as of 2025 January 12. Given that GLIMPSE adopts a subpixel dithering pattern to minimize overheads, any impurities in the flats can be coadded in the dithering procedure and thus might appear as real sources in each primary dither position. Compared to the publicly available flats provided by the Space Telescope Science Institute available as of 2025 mid-January (first released in 2023 September), our flats result in substantially improved depth in the final mosaics. We recover an ≈0.3–0.5 mag improvement in depth across all long-wavelength (LW) bands and an ≈0.1–0.2 mag improvement across the SW bands. Finally, given the depth of the GLIMPSE campaign and the presence of bright cluster galaxies (bCGs), we model and subtract the background on an amplifier basis with SEP (K. Barbary 2016), while manually masking out the bright regions. This further improves the image background and allows the mosaics to reach ~0.3 mag deeper across all filters.

Our final data set achieves 5σ aperture-corrected nominal depths of 30.8–30.9 mag across all broad bands in $D = 0''2$ apertures.

We also process and incorporate the deep HST ACS and WFC3 mosaics from the HFF (J. M. Lotz et al. 2017) and BUFFALO (C. L. Steinhardt et al. 2020) programs. The HST images are based on Gaia-aligned mosaics from the CHArGE archive (V. Kokorev et al. 2022), which are hosted on the Dawn JWST Archive (F. Valentino et al. 2023). We drizzle our final mosaics onto a $0''02 \text{ pixel}^{-1}$ grid for the JWST SW filters and $0''04 \text{ pixel}^{-1}$ for JWST LW and HST.

## 2.3. Cluster Light Contamination

Careful handling of the contamination light from the bCGs and intracluster light (ICL) is necessary for our science objectives. Improper treatment of the ICL could potentially leads to inaccurate colors, which in return affect the derived photometric redshifts, negatively impacting the high-z selection in the proximity of the cluster. To model and subtract the bCG and ICL light, we follow the robust methods described in L. Ferrarese et al. (2006), H. V. Shipley et al. (2018), and J. R. Weaver et al. (2024) for the Hubble Frontier Fields-DeepSpace (HFFDS) and UNCOVER/Mega Science photometric catalogs (K. A. Suess et al. 2024). After both the bCGs and the ICL have been modeled and subtracted, we perform an additional local background subtraction pass in the affected areas, largely following the same methodology as in Section 2.2.

## 2.4. Source Extraction

We construct empirical PSFs by stacking the available stars in the field, and we match all available bCG-subtracted JWST and HST images to our lowest-resolution PSF—F480M. The sources are detected using SExTRACTOR (E. Bertin & S. Arnouts 1996), which we run in a single mode on an inverse-variance-weighted combination of the native PSF F277W, F356W, and F444W images. We then use PHOTUTILS (L. Bradley et al. 2020) to measure flux densities in circular apertures with varying diameters from $D = 0''1$ to $1''2$. The photometric errors are determined separately for each object and filter, taking into account both the aperture size and the depth variation across the image. In the vicinity of each source, we place 2000 random apertures in source-free parts of the image (as determined by the segmentation map). The standard deviation of the flux density within empty apertures plus the Poisson noise are then used as the final uncertainty. This method is generally preferable over just using the weight/error maps, as it better accounts for effects of correlated noise (see, e.g., R. Endsley et al. 2023; J. R. Weaver et al. 2024). The aperture corrections are calculated by assuming a point-source profile. We use the empirical F480M symmetric PSF curve of growth to determine the fraction of the total flux that falls outside of each specified aperture size and then use that as our correction factor.

In order to select high-fidelity high-z candidates, we flag all sources that fall on the edge of the mosaic, intersect with diffraction spikes from bright stars, or are close (<2''0) to the modeled and subtracted bCGs. As galaxies at high $z$ are generally small, we will only use the total fluxes within a $D = 0''2$ aperture for the remainder of our work.

# 3. Data Analysis

Our goal is to find galaxies above $z > 16$. Previous experience has demonstrated that a single selection technique, however, is insufficient to provide robust $z > 16$ candidates and might result in low-z contaminants (see, e.g., R. P. Naidu et al. 2022). To mitigate that, we combine our selection to use the Lyman break (Section 3.1) and photometric redshift test (Section 3.2), which stringently tests for low-z contamination (Section 3.3) and uses the size estimates (Section 3.6) to further test their robustness.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** High-z color–color selection. Our final high-z sample, which we also present in Figure 3, is shown with red diamonds. The color–color Lyman break selection is indicated with a black dashed line. Solid lines show the tracks followed by starburst galaxies at $z \sim 16$ (J. Chevallard & S. Charlot 2016) with varying levels of attenuation $A_V = 0-0.5$ (shown with blue to green). We highlight low-redshift quiescent galaxies (M. Polletta et al. 2007), whose Balmer break could mimic the Lyman break at higher redshifts, as dashed lines. Extreme dusty starbursts (R. P. Naidu et al. 2022; P. Arrabal Haro et al. 2023; J. A. Zavala et al. 2023) at $z \sim 4.5-5.5$ (orange squares) look quite similar to a high-z galaxy with this selection alone (see Section 3.3 for why these are disfavored). Finally, cool stars and brown dwarfs (G. Chabrier et al. 2000; F. Allard et al. 2001) are shown in purple.</figcaption></figure>

### 3.1. Lyman Break Selection

The first selection applied to the photometric catalog is based on identifying dropouts in the rest-frame UV, using the $Ly\alpha$ break technique. We define a selection window in color–color space to isolate dropouts between the F277W and F200W filters, while simultaneously excluding objects with red continua in the rest-frame UV between the F277W and F356W filters. These color criteria are determined from running synthetic photometry on a set of galaxy templates generated using BEAGLE (J. Chevallard & S. Charlot 2016). We used starburst templates in the redshift range $15 < z < 20$, incorporating intergalactic medium attenuation following A. K. Inoue et al. (2014), with varying levels of attenuation between $A_V = 0$ and $A_V = 0.5$, based on the SMC extinction law (K. D. Gordon et al. 2003). We also consider low-redshift quiescent galaxies that could mimic a Lyman break, using spectral energy distributions (SEDs) from the SWIRE template library (M. Polletta et al. 2007). In addition, we include cool stars and brown dwarfs, using stellar templates from G. Chabrier et al. (2000) and F. Allard et al. (2001). The color–color tracks for all these simulated sources are shown in Figure 2, along with the adopted selection window, which is defined by the following criteria:

$$ 
\begin{aligned}
M_{200} - M_{277} &> 1.0 \\
M_{200} - M_{277} &> 1.2 + 2.0(M_{277} - M_{356}) \\
M_{277} - M_{356} &< 0.5.
\end{aligned}
$$

In addition to the color criteria, we require that sources are well detected (at the $>3\sigma$ level in all three of the broadband

LW filters) and have an $S/N$ of $\gtrsim 4$ in at least one band, while remaining undetected (at the $<2\sigma$ level) in bands blueward of the Lyman break (F090W, F115W, and F150W). This first selection leads to a sample of $\sim 38 z \gtrsim 16$ galaxy candidates.

### 3.2. Photometric Redshift

To calculate photometric redshifts ($z_{phot}$) for all objects in the GLIMPSE catalog, we use the PYTHON version of EAZY (G. B. Brammer et al. 2008). We choose the BLUE_SFHZ_13 model subset$^{20}$ that contains redshift-dependent SFHs and dust attenuation values. More specifically, the linear combinations of lognormal SFHs included in the template set are not allowed to exceed redshifts that start earlier than the age of the Universe (for more details see M. R. Blanton & S. Roweis 2007). These models are further complemented by a blue galaxy template, derived from a JWST spectrum of a $z = 8.50$ galaxy with extreme line equivalent widths (A. C. Carnall et al. 2023).

Our work is concerned only with the F200W dropouts, for which the HST data are simply too shallow to be useful. Therefore, in our analysis we only focus on the GLIMPSE JWST bands. We fit the aperture-corrected $D = 0.2$ JWST flux densities, including the upper limits, using the $0.01 < z < 30$ redshift grid. The uncertainties on the photometric redshifts are computed from the 16th and 84th percentiles of the redshift probability distributions—$p(z)$. The best-fit EAZY SEDs are only used to validate the color–color high-z selection.

In order to select $z > 16$ galaxy candidates with EAZY, we then require that the sources must be detected ($S/N > 3$) in at least three bands, similarly to our color–color selection, and have a well-constrained $p(z)$ (FWHM $< 2.5$) without a statistically significant secondary redshift solution. Cross-referencing our EAZY high-$z$ sample with the color–color-selected F200W dropout from Section 3.1 results in eight candidates at $z \simeq 15.6-19.5$.

In addition to EAZY, we use BEAGLE (J. Chevallard & S. Charlot 2016) and the methods outlined in R. Endsley et al. (2024) to independently derive photometric redshifts. Briefly, we use the BEAGLE models constructed from updated G. Bruzual & S. Charlot (2003) stellar templates that use the PARSEC isochrones (A. Bressan et al. 2012; Y. Chen et al. 2015). These models have been passed through Cloudy (G. J. Ferland et al. 2017) to self-consistently produce nebular (both line and continuum) emission (J. Gutkin et al. 2016). We explore different parametric SFHs (including constant, bursty, and a two-component SFH), stellar ages (between 1 Myr and 30 Gyr), metallicities (0.0063–0.5 $Z_{\odot}$), and attenuation laws (we settle on the SMC law; Y. C. Pei 1992). While the two codes agree exceptionally well (Figure 3) for the majority of our high-$z$ candidates, we further remove two sources where the BEAGLE fit shows a prominent secondary peak solution, leaving us with six candidates.

### 3.3. Low-z Confusion

In the early days of JWST, initial discoveries saw a slew of $z > 16$ galaxy candidates, for which valid concerns have been raised in the literature regarding the ability of photometric redshifts to pick out true high-$z$ candidates (P. Arrabal Haro et al. 2023; C. T. Donnan et al. 2023; J. A. Zavala et al. 2023). These works have shown that dusty extremely star-forming galaxies at $z \sim 5$ can masquerade as high-$z$ objects and bias our

<footnote>$^{20}$ https://github.com/gbrammer/eazy-photoz/tree/master/templates/sfhz</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Most secure *z* ≳ 16 candidates from GLIMPSE. For each source we show 1''5 cutouts in all broad- and medium-band NIRCam filters, plus a detection (F277W +F356W+F444W) LW stack. On each cutout we overlay the extraction aperture (*D* = 0''2) in green. On the SED plot, we show the best-fit high-*z* vs. low-*z* EAZY models (maroon and blue, respectively), as well as the same for BEAGLE (dashed orange and green). The colored circles represent the integrated flux density for each solution. Detections are shown as black squares; 3σ upper limits are shown as downward-pointing arrows. We explicitly highlight the difference between a low-*z* interloper and high-*z* EAZY solution with a ΔBIC statistic in each panel, as discussed in Section 3.3. In the inset panel we show the *p(z)* for both EAZY and BEAGLE fits with a free redshift solution. Finally, we display the probability of these candidates being at *z* > 14 as derived from multiplying EAZY and BEAGLE redshift distribution.</figcaption></figure>
view of high-redshift galaxies. For example, the so-called Schrödinger galaxy was initially photometrically identified at $z > 16$ (R. P. Naidu et al. 2022; C. T. Donnan et al. 2023) and then had its redshift revised to $z \sim 5$, via extreme emission lines found in NIRSpec spectra, as an interloper (P. Arrabal Haro et al. 2023; J. A. Zavala et al. 2023). Despite this misidentification, these outliers provide invaluable insight into the potential pitfalls of our methodology. Subsequently, identification methodologies have been improved, and the techniques used to process (e.g., various pipeline improvements and calibrations) and analyze JWST data have significantly matured, such that the vast majority of $\sim 95$ photometrically identified $z > 10$ galaxies have been spectroscopically confirmed at high redshift (P. Arrabal Haro et al. 2023; K. N. Hainline et al. 2024; Y. Harikane et al. 2025). At this moment, these initial misidentifications provide a road map to rigorously scrutinize the $z > 16$ candidates; however, at the moment no galaxies with $z > 15$ have been confirmed spectroscopically.

To address the potential contamination from low-$z$ interlopers, we refit all nine of our high-$z$ candidates with the EAZY BLUE_SFHZ_13 model suite, plus an additional dusty starburst template (R. P. Naidu et al. 2022; P. Arrabal Haro et al. 2023), to mimic the low-$z$ interlopers. The redshift grid in this case is forced to be $z < 6$, to match the redshift range of dusty interlopers. Deciding which model is better relies on the Bayesian information criterion (BIC) test (G. Schwarz 1978). This test computes which template is the most likely fit to the observed data, while penalizing models that have too many free parameters. In our case a high-$z$ fit is statistically preferred over forced low-$z$ ($z < 6$) solutions with the additional Schrödinger template, when a BIC difference reaches 3 or more (using the criteria defined in H. Jeffreys 1961). We find that a majority ($5/6$) of our high-$z$ objects can be fit equally well with a low-$z$ dusty starburst template ($\Delta$BIC $< 3$, corresponding to $\gtrsim 3\sigma$ significance). We are now left with five $z \sim 16$ galaxy candidates. While these appear to be robustly detected ($\sim 3\sigma$) in three bands, we further limit our final sample to only include galaxy candidates that reach at least $5\sigma$ in one or more bands. Doing this further removes three objects that are located on the edges of our mosaics, where we do not have the full coverage of all six primary dithers (see Figure 1) and thus do not achieve the full depth. We further verify that our remaining objects are not simply hot pixels in the LW bands, by separately examining all six primary dithers in the F277W band, as well as the LW stack.

This final quality cut ensures that only the best high-$z$ candidates end up in our final sample. This leaves us with two objects, which we show in Figure 3. It is worth pointing out that the secondary $z < 6$ solution with a Schrödinger template is choosing exactly the redshift range ($z \sim 4.9$) where emission lines can confuse the redshift-fitting codes (R. P. Naidu et al. 2022; P. Arrabal Haro et al. 2023); however, given our chosen $\Delta$BIC threshold, high-$z$ EAZY solutions are still preferred. To further test this, we carry out an additional BEAGLE low-$z$ fit where we very finely sample the $z = 4.8-5.1$ range in an attempt to find the highly specific redshift solution where high emission line intensities can mimic a break. Our results find that the finer spacing grid does not change the conclusion in a statistically significant way. Finally, we note that the morphology of our objects (Figure 3) does not change from band to band when moving to redder filters, as it would for dusty galaxies.

### 3.4. Transient Contamination

Low-$z$ confusion could also arise from transient events, as was noted in C. DeCoursey et al. (2025), who showed that certain supernovae can mimic SEDs of $z > 16$ galaxies. Generally, multiepoch observations with long-enough baselines can rule out these contaminants; however, GLIMPSE data with only $\sim 24$ hr between LW observations are not suitable to detect such changes. Despite that, as C. DeCoursey et al. (2025) point out, deeper observations in SW channels can aid in distinguishing a supernova from a Lyman break, as the drop-off of the former is more gradual. Given the sharp breaks, aided by deep F200W data, and extended morphologies of our candidates (Figure 3), contamination by supernovae seems unlikely.

### 3.5. Gravitational Magnification

To take into account the effects from gravitational lensing, we use a new strong-lensing (SL) model of AS1063 constructed for GLIMPSE with the updated version of the A. Zitrin et al. (2015) parametric method that was recently used for several JWST SL clusters (e.g., M. Pascale et al. 2022; L. J. Furtak et al. 2023; A. K. Meena et al. 2023). We model the cluster with two smooth DM halos parameterized as pseudoisothermal elliptical mass distributions (A. Kassiola & I. Kovner 1993): one centered on the bCG, and the other on a group of galaxies in the northeast of the cluster (e.g., P. Bergamini et al. 2019; B. Beauchesne et al. 2024). In addition, we model 303 cluster member galaxies as dual pseudoisothermal ellipsoids (Á. Elásdóttir et al. 2007). The model is constrained with 80 multiple images of 30 sources, 26 of which have spectroscopic redshifts (I. Balestra et al. 2013; A. Monna et al. 2014; J. Richard et al. 2021; B. Beauchesne et al. 2024; M. W. Topping et al. 2024). The model achieves an average image reproduction error of $\Delta_{rms} = 0.47$. We refer the reader to L. Furtak et al. (2025, in preparation) for more details on the lens model. A previous (HST-based) version of this model was recently used in M. W. Topping et al. (2024).

Magnifications are computed analytically at the position of each galaxy candidate and adopting its photometric redshift. The magnifications and their uncertainties are listed in Table 1.

### 3.6. Size Measurements

The cutouts we show in Figure 3 imply that our sources are quite compact, yet they appear to be resolved. To measure the effective radii, we model each source with GALFIT (C. Y. Peng et al. 2002, 2010), using a Sérsic (J. L. Sérsic 1963) profile where the source position, brightness, effective radius, Sérsic index, and axis ratio are allowed to vary. When fitting, we take into account the effects of the PSF, which we measure empirically from the bright stars in the field. We perform this procedure on our brightest band, F277W, to ensure that optimal S/N per pixel is achieved to accommodate robust size measurements. Additionally, in the case of source 70,464, we simultaneously model both the high-$z$ object and the low-$z$ ($z_{phot} \sim 3$) neighboring galaxy. We find that the on-the-sky sizes of our sources are in the range of $R_{eff} \sim 0.08 - 0.10$. A source can be considered to be resolved when its effective radius is larger than the empirical PSF half-width at half-maximum (HWHM). Since the HWHM of F277W PSF is $\sim 0.046$, we can consider all our sources to be resolved.

Taking the redshift and gravitational magnification into account, we convert our angular sizes to physical effective
<table><caption>**Table 1**
Properties of the High-z Sample<sup>a</sup></caption><thead><tr><th>Parameter</th><th>70467</th><th>72839</th></tr></thead><tbody><tr><td>R.A. (deg)</td><td>342.2386</td><td>342.2319</td></tr><tr><td>Decl. (deg)</td><td>-44.5347</td><td>-44.5565</td></tr><tr><td>z<sub>phot</sub></td><td>16.4 ± 1.8</td><td>15.8 ± 0.8</td></tr><tr><td>Δχ<sup>2</sup></td><td>1.0</td><td>1.2</td></tr><tr><td>ΔBIC<sup>b</sup></td><td>4.1</td><td>3.3</td></tr><tr><td>μ</td><td>1.55 ± 0.04</td><td>1.35 ± 0.02</td></tr><tr><td>M<sub>UV</sub> (ABmag)</td><td>-17.0 ± 0.2</td><td>-17.2 ± 0.2</td></tr><tr><td>SFR (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>0.4 ± 0.1</td><td>0.5 ± 0.1</td></tr><tr><td>β</td><td>-2.91 ± 0.20</td><td>-2.82 ± 0.33</td></tr><tr><td>R<sub>eff</sub> (pc)</td><td>212 ± 23</td><td>199 ± 71</td></tr></tbody></table>

**Notes.**

<sup>a</sup> All values are corrected for gravitational magnification.

<sup>b</sup> BIC(low-z) - BIC(high-z).

radii. These are roughly similar for both objects with $R_{\text{eff}} \sim 200$ pc. We list the delensed sizes in Table 1.

# 4. Results

## 4.1. Stellar Population Properties

After completing the multiple stages of selecting our final $z > 16$ sample, we now compute a range of relevant physical parameters. The uniqueness of our sample lies in the unprecedented depth of the GLIMPSE observations, which we further push to their limit with gravitational lensing in order to obtain these candidates. On average, our sources have three individual band detections, with the rest of the JWST photometry being upper limits. While this is perfectly adequate to constrain the redshift, deriving any physical parameters, especially stellar mass, from SED fitting codes would be simply unreasonable. As such, we will only focus on the observable stellar population parameters that can be derived from the rest-UV photometry alone.

We estimate the absolute UV magnitude for each galaxy from the observed F277W band, which samples $\lambda_{\text{rest}} \sim 1500$ Å at this redshift. We find that our targets cover a very narrow range of derived $M_{\text{UV}} = -17.1_{-0.12}^{+0.10}$, after accounting for lensing magnification. We derive the UV slope $\beta$ for each object by assuming $f_{\lambda} \sim \lambda^{\beta}$ and fitting it to all the photometric points that fall within the $\lambda_{\text{rest}}$ in the 1260–2500 Å range, which effectively traces the observed F277W–F356W color. We find that our $\beta$ values span a very narrow range from $-2.9$ to $-2.8$, with a median of $-2.87 \pm 0.15$. These are shown in Figure 4, alongside the magnification-corrected UV brightness. Finally, we derive the SFR$_{\text{UV}}$ directly from our delensed $M_{\text{UV}}$, by following the relation from R. C. Kennicutt & N. J. Evans (2012). We find that the SFR$_{\text{UV}}$ ranges from 0.4 to 0.6 $M_{\odot}$ yr$^{-1}$. Since our derived $\beta$ values imply negligible dust reddening, we do not apply a correction for dust when deriving the star formation rate (SFR). All the delensed physical parameters are listed in Table 1.

## 4.2. Number Density

We compute the number density of $16 < z < 20$ sources in GLIMPSE by utilizing a $1/V_{\text{max}}$ method (M. Schmidt 1968), where $V_{\text{max}}$ corresponds to the maximum volume a galaxy could occupy and still be detected in the appropriate filter for our redshift range. Due to lensing, the effective area, and therefore volume, covered by our observations is smaller than it

<figure><img src="image_6.png" /><figcaption>**Figure 4.** UV continuum slope and delensed luminosity. The black dashed line shows the best-fitting $\beta - M_{\text{UV}}$ relation from F. Cullen et al. (2023) for star-forming galaxies at $z = 8-16$. The gray shaded area corresponds to the 68% confidence interval. Our secure ($>5\sigma$) sources at $z \sim 16$, shown with red diamonds, appear to be consistent with the predicted trend within $1\sigma$, although systematically shifted to bluer slopes.</figcaption></figure>

would be in a blank field. In addition, an accurate derivation of the UV number density requires a robust estimate of the completeness of our survey, for a given selection function and depth. To account for the former, we use our magnification maps together with the $p(z)$ for our sources (to account for the redshift range) and, including the uncertainties on both, derive a total survey volume of $17.055 \pm 71 \text{ Mpc}^3$.

The second step to derive the effective volume consists of computing the survey completeness through the lensing cluster. To do that, we use the same approach adopted in H. Atek et al. (2018) and I. Chemerynska et al. (2024). Briefly, the procedure includes generating a large set of mock galaxies that are distributed directly in the source plane, which in turn was generated from our lensing model. The properties of the simulated objects span the redshift, size, and luminosity range of our objects. More details regarding this procedure will be presented in I. Chemerynska et al. (in preparation). Combining this with our derived volume, we find an effective survey volume $V_{\text{max}} = 5882.7 \pm 1538.9 \text{ Mpc}^3$. The final uncertainty on the volume is derived from the combined errors on the lensing model (to be presented in L. Furtak et al. 2025, in preparation), the redshift probability distribution of the high-z candidates, and the completeness simulation itself.

All of our objects cover nearly the same range of $M_{\text{UV}}$, so to derive the final number density we will just assume a single bin with a width of 1 mag—$M_{\text{UV}} = -17.0 \pm 0.5$. The final number density for our sample in the $16 < z < 20$ range is therefore equal to $\log_{10}(\phi/[\text{Mpc}^{-3}\text{mag}^{-1}]) = -3.47_{-0.10}^{+0.13}$.

# 5. Discussion

## 5.1. Stellar Populations at High Redshift

In this section, we reflect on the implications of our derived parameters if the selected sample truly resides at $z \sim 16$.

First, we compare our calculated $R_{\text{eff}}$ with empirical results for spectroscopically confirmed galaxies at $z \sim 13-14$. Both
S. Carniani et al. (2024b) and E. Curtis-Lake et al. (2023) report UV sizes in the range of ~100–300 pc, consistent with our findings. While high-redshift size predictions from simulations are limited, the TNG50 simulations have accurately reproduced galaxy morphologies across a wide redshift range (see, e.g., S. Tacchella et al. 2019). The latest high-z size evolution analyses in L. Costantin et al. (2023) and T. Morishita & M. Stiavelli (2023) suggest a scaling relation of ~8.66 × (1+z)⁻¹.¹⁵ kpc. Extrapolating this to z ~ 16 predicts R<sub>eff</sub> ~ 250–300 pc, aligning with our measurements within 1σ. Generally, the relatively large sizes of our objects may indicate that UV light from these galaxies arises from extended stellar populations, as observed in S. Carniani et al. (2024b). Furthermore, the spatially resolved nature of our sources, combined with the absence of a PSF-like, centrally concentrated component, suggests that the UV emission is not dominated by unobscured active galactic nuclei (AGN), unlike other high-z galaxies (Y. Harikane et al. 2023; R. Maiolino et al. 2024).

Using our derived SFR<sub>UV</sub> and R<sub>eff</sub>, we calculate the UV-based SFR surface density, finding a median value of Σ<sub>SFR,UV</sub> = 1.14<sup>+0.20</sup><sub>-0.18</sub> M<sub>☉</sub> yr<sup>-1</sup> kpc<sup>-2</sup>. This is 5–10 times lower than values reported for spectroscopically confirmed objects at z ~ 14 by S. Carniani et al. (2024b) and is similarly lower when extrapolating the trend in A. Calabrò et al. (2024). However, higher Σ<sub>SFR,UV</sub> found at z < 15 likely reflect observational biases and small number statistics at high z, with the bulk of the population likely containing less intense star formation. Moreover, if the accretion rate changes exponentially with redshift for galaxies such as ours, even small changes in redshift will likely lead to significantly more accretion. We will explore this further in Section 5.3.

The UV slopes in our sample also do not exhibit unusual characteristics. We find a median $\beta$ of -2.87, consistent with a largely dust-free stellar plus nebular continuum. This is reasonable, as our objects are observed just ~200 Myr post-big bang, with limited time for substantial dust production. Within the uncertainties of our $\beta$ values, our sources align with the $\beta$-$M_{UV}$ luminosity relation from the surveys used in F. Cullen et al. (2023), which are based on recent JWST data for z = 8–16 star-forming galaxies. Finally, we do not observe extremely negative beta slopes that have been postulated for metal-free Population III stars (see, e.g., J. Jaacks et al. 2018).

## 5.2. Where Are the Bright High-z Galaxies?

The GLIMPSE imaging data are contained within a single NIRCam pointing, with the effective area further limited by lensing effects. Despite the relatively small volume surveyed, we identify two robust high-z candidates, all exhibiting characteristics consistent with a star-forming population at this epoch. We now discuss the number density derived from our sample, comparing it with theoretical predictions from both pre- and post-JWST launch models (see Figure 5).

Regardless of the simulation type—whether semianalytic (P. Dayal et al. 2014; C. A. Mason et al. 2015; P. Behroozi et al. 2019) or hydrodynamical (J. Rosdahl et al. 2022; S. M. Wilkins et al. 2024)—we find that our derived number density at 16 < z < 20 and $M_{UV} = -17.0 \pm 0.5$ significantly exceeds the pre-JWST theoretical expectations, with a discrepancy of over $3\sigma-5\sigma$. This mirrors an observed discrepancy in the UV luminosity function (UVLF) suggested by spectroscopically confirmed z ~ 12–14 objects (C. T. Donnan et al. 2023; S. L. Finkelstein et al. 2024; S. Carniani et al. 2024b;

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Theoretical comoving number densities of $z > 16$ galaxies in the literature. Our completeness-corrected densities are shown in red. Black lines are median values; shaded areas show 68th percentiles. Values have been homogenized in terms of redshift interval (16 < z < 20) and delensed. While high, our predictions generally align with the upper edge of the theoretical predictions (P. Dayal et al. 2014, 2019; C. A. Mason et al. 2015; P. Behroozi et al. 2019; V. Mauerofer & P. Dayal 2023; J. B. Muñoz et al. 2023; R. Feldmann et al. 2025; Z. Li et al. 2024).</figcaption></figure>

C. M. Casey et al. 2024; I. Chemerynska et al. 2024; Y. Harikane et al. 2025; B. Robertson et al. 2024). In contrast, post-JWST calibrated simulations, including semianalytic models (V. Mauerofer & P. Dayal 2023; J. B. Muñoz et al. 2023), analytic predictions that include the feedback-free burst scenario (Z. Li et al. 2024), and FIREbox hydrodynamical simulations (R. Feldmann et al. 2025), show better (generally within $1\sigma-2\sigma$) agreement.

We find the best agreement with the FIREbox simulations presented in R. Feldmann et al. (2025). This might suggest that the abundance of UV-luminous galaxies at $z \gtrsim 11$ (P. Arrabal Haro et al. 2023; E. Curtis-Lake et al. 2023; R. Maiolino et al. 2024; S. Carniani et al. 2024b; M. Castellano et al. 2024) could be due to a relatively constant star formation efficiency (SFE) that is largely independent of $M_{halo}$ (I. Chemerynska 2025, in preparation). We will explore the connection of our sample to these lower-z UV-bright galaxies in the next section. Encouragingly, our sample, while showing somewhat higher number densities at this absolute magnitude, remains broadly consistent with many post-JWST theoretical predictions.

Over the past 2 yr, JWST has identified dozens of bright $M_{UV} < -18$ high-z galaxy candidates at $z = 10-14$, challenging pre-JWST expectations (N. J. Adams et al. 2023; H. Atek et al. 2023; D. Austin et al. 2023; E. Curtis-Lake et al. 2023; C. T. Donnan et al. 2023; S. L. Finkelstein et al. 2024; G. C. K. Leung et al. 2023; P. G. Pérez-González et al. 2023; S. Carniani et al. 2024b; M. Castellano et al. 2024; I. Chemerynska et al. 2024; D. J. McLeod et al. 2024; B. Robertson et al. 2024). Surprisingly, despite the large areas and depths of various surveys, none have detected bright
galaxies beyond $z > 16$. Given the detection limits of ~30–30.5 mag in the F200W/F277W filters (used for Lyman break dropout selection at $z > 15$), bright galaxies with $M_{UV} < -19$ would likely have been identified by these surveys if they indeed existed at these redshifts. Selecting these galaxies at $z > 16$ requires deep F200W (dropout band) and F277W photometry, as the break has to be unambiguous ($>1$ mag). To compute the expected number densities from various extragalactic surveys, based on nondetections (N. Gehrels 1986), we do the following. Using the quoted F200W $2\sigma$ depths for CEERS (M. B. Bagley et al. 2023; S. L. Finkelstein et al. 2024), PRIMER (both COSMOS and UDS), JADES Origins Field (D. J. Eisenstein et al. 2023; B. Robertson et al. 2024), and NGDEEP (G. C. K. Leung et al. 2023; M. B. Bagley et al. 2024), we assume $M_{200} - M_{277} \ge 1$ mag and then compute the expected $M_{UV}$ and the corresponding number density (from $V_{max}$) for each survey using the reported survey area. These are shown in Figure 6 as gray upper limits, alongside the GLIMPSE number density at $M_{UV} = -17$ and upper limits at $M_{UV} = -19$ and $-18$. The discrepancy between the limits at $M_{UV} < -18$ and our observed number densities at $M_{UV} = -17$ is quite jarring. The upper limits reach ~$10^{-6}$ Mpc$^{-3}$ mag$^{-1}$ (for CEERS plus PRIMER) at $M_{UV} \sim -19$, in a stark contrast to GLIMPSE, where we predict that there are $\gtrsim 2000 \times$ more galaxies per volume at $M_{UV} = -17$. If $M_{UV} \lesssim -18$ galaxies existed at this epoch, they surely should have been detected by now.

So where are they? One possibility is that because of cosmic variance, the current deep fields have missed some of the bright $z > 15$ galaxies. After all, the depths of GLIMPSE are comparable to JADES and NGDEEP; as such, we would expect these surveys to find a $z > 15$ source. In this scenario, wider-field surveys that reach $M_{UV} \sim -20$ mag would discover bright $z > 15$ galaxies. For instance, objects brighter than $M_{UV} \sim -18$ are too rare to appear in the FIREbox-simulated volume at this redshift (see Figure 6).

Furthermore, this difference can be potentially due to the bright end of the UVLF at $z > 16$ lying exactly at the GLIMPSE detection limit of $M_{UV} \sim -17$, explaining why we are seeing these objects for the first time (though see a candidate in NGDEEP found by both G. C. K. Leung et al. 2023 and D. Austin et al. 2023 at $z_{phot} \sim 15.6$, with $M_{UV}$ of $\sim -19.2$, which implies a $\log_{10}(\phi/[\text{Mpc}^{-3}\text{mag}^{-1}]) \sim -4.65$). We would also like to further highlight several relatively bright ($M_{UV} \sim -19$) $z > 16$ photometric candidates identified in K. N. Hainline et al. (2024) and C. J. Conselice et al. (2024), which suggest comparable number densities to those found in G. C. K. Leung et al. (2023) and D. Austin et al. (2023). However, these objects exhibit weaker Lyman breaks than the candidates in our study. Before any of these objects are confirmed spectroscopically, the fate of the UVLF at high $z$ would remain elusive.

The steepness of the UVLF implied by the GLIMPSE detections against the brighter limits can be fit if the UVLF closely tracks the halo mass function, as we show in Figure 6 (translating halo mass to SFR with a toy model given by the cosmic baryon fraction, zero dust, and 30% efficiency of gas into stars over a fixed 100 Myr timescale). We note that the most plausible way to reproduce the observed number density is with a constant and relatively high SFE.

The lensing magnification present in AS1063, as opposed to deep blank fields such as JADES and NGDEEP, could also play a role in this discrepancy. We speculate that even the modest lensing magnification of $\mu \sim 1.5$, which will double the

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Observed abundance of bright galaxies at $z > 16$. We report the approximate magnitude limits for a number of deep extragalactic surveys corresponding to the $2\sigma$ depth of the dropout filter—F200W, which is required for the Lyman break identification. In varying shades of gray, we show CEERS (M. B. Bagley et al. 2023; S. L. Finkelstein et al. 2024) and PRIMER (both COSMOS and UDS; C. T. Donnan et al. 2023) combined, JADES (Origins Field; D. J. Eisenstein et al. 2023), and NGDEEP (G. C. K. Leung et al. 2023; M. B. Bagley et al. 2024). The dashed line is a scaled version of the halo mass function that assumes a 30% SFE of the gas, no dust, and a continuous SFH. The dashed line is simply a toy model to illustrate how the UVLF at these redshifts could steeply decline.</figcaption></figure>

exposure time on source (adding ~0.5 mag in depth), will provide the additional sensitivity needed to uncover these early galaxies. There simply was no feasible way to see these relatively rare objects prior to GLIMPSE.

### 5.3. Abundance of Bright Galaxies at $z \sim 12-14$

Our observations indicate that one possible explanation for the lack of $M_{UV} \lesssim -18$ galaxies at $z \sim 16$ is that brighter galaxies are extremely rare at high $z$, so only “faint” galaxies would be abundant enough to be found in the volume probed so far at this epoch. Therefore, one can imagine that the overabundance of bright galaxies observed at $z = 12-14$ has likely evolved from earlier (f fainter) galaxies, such as the ones presented here, which are surprisingly abundant for their relatively faint UV luminosities. In return, this suggests that the discrepancy between theoretical predictions and observations may begin as early as $z = 18$, potentially during the epoch when the first galaxies formed.

To test that, we would like to see how different SFHs starting at $z \sim 16$ reproduce the observed bright galaxies at $z < 14$. To do this, we consider two cases. First, we consider a constant SFH that begins forming stars at the median redshift of our sample ($z \sim 16.1$) at a rate of $\sim 1 \text{ M}_\odot \text{yr}^{-1}$ (to match our derived SFRs). As the second case, we would consider the maximum growth rate allowable for these galaxies by forming stars at 100% of the accretion rate set by the extended Press–Schechter approximation in a ΛCDM Universe (see, e.g., Equation (7) in A. Dekel & M. R. Krumholz 2013, with a 0.3 dex scatter as in K. Ren et al. 2018 and J. Mirocha et al. 2021). This in turn creates an exponentially increasing SFH that puts a maximum limit on the luminosity growth experienced by these galaxies observed at these redshifts. Since we want to show the absolute limit possible, given our observations, we start to evolve this SFH at the 84th percentile ($1\sigma$) of the median redshift of our
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The abundance of bright sources at cosmic dawn. We trace the possible evolutionary paths of our objects from $z \sim 16$, while factoring in the redshift uncertainty, to some of the brightest spectroscopically confirmed high-$z$ objects at $z \simeq 11–14$ (P. Arrabal Haro et al. 2023; A. J. Bunker et al. 2023; E. Curtis-Lake et al. 2023; S. Carniani et al. 2024b; M. Castellano et al. 2024). Shown are the tracks (with $1\sigma$ uncertainty represented by shaded regions) of constant star formation at $\sim 1 M_{\odot} \text{yr}^{-1}$ (purple) and a $\Lambda$CDM-limited maximal possible accretion scenario (pink line; A. Dekel & M. R. Krumholz 2013), where all gas is converted to stars. The overabundance of bright galaxies observed at high $z$ so far is fully consistent with our observations. We only show the $M_{UV}$ uncertainties that are larger than the marker size.</figcaption></figure>

sample, $z \sim 18.0$, instead of simply the median. Both scenarios begin at $M_{UV} \sim -17$. To evolve these SFHs, we have used the tool kit within BAGPIPES (A. C. Carnall et al. 2019), which outputs $M_{UV}$ in a predefined grid of look-back time. We show the results of our modeling, alongside the bright spectroscopically confirmed galaxies, in Figure 7.

We find that if our galaxies followed a constant SFH and continued forming stars at their current rate, they can easily reproduce the objects observed in E. Curtis-Lake et al. (2023) and the fainter one of the two galaxies presented in S. Carniani et al. (2024b). Reproducing some of the brightest ($M_{UV} < -20$) known galaxies at $z \ge 10$ would, however, require very (nearly 100%) efficient star formation episodes. Following the maximal exponential accretion scenario, we find that galaxies in our sample can, within $1\sigma$, become as bright as JADES-GS-z14-0 ($M_{UV} = -20.81$ at $z = 14.17933$; S. Carniani et al. 2024b, 2024a; S. Schouws et al. 2024) or GNz11 (although the stellar contribution to the $M_{UV}$ could be lower, due to the potential AGN; P. A. Oesch et al. 2016; A. J. Bunker et al. 2023; R. Maiolino et al. 2024), without a need to invoke any exotic cosmological framework. Curiously, the SFH for JADES-GS-z14-0 presented in S. Carniani et al. predicts $\text{SFR} \lesssim 1 M_{\odot} \text{yr}^{-1}$ at $z \sim 16–17$, which matches our sample if it were a potential progenitor. This further reemphasizes that these $z \sim 16$ candidates could plausibly grow to become the extremely bright $z \sim 14$ galaxies without extraordinarily new cosmologies.

## 6. Conclusions

In this study, we leverage ultradeep GLIMPSE NIRCam imaging covering the AS1063 lensing cluster to identify a unique sample of high-redshift ($z > 16$) galaxy candidates. Using a combination of observed colors alongside an extensive

suite of SED modeling routines, we select two robust candidates at $z \simeq 15.8–16.4$, which all fit into our single NIRCam pointing.

We examine the stellar population properties of our sample, such as their delensed, absolute UV brightness, $\beta$ slopes, and sizes, finding that these high-$z$ galaxies fit very well into the picture of galaxy evolution being established at $z \sim 10–14$. Crucially, our findings reveal a significant excess in number density of galaxies at $M_{UV} = -17$ of $\log_{10}(\phi/[\text{Mpc}^{-3}\text{mag}^{-1}]) = -3.47_{-0.10}^{+0.13}$, compared to theoretical predictions from both semianalytic and hydrodynamical simulations calibrated before the JWST launch, mirroring similar discussions for bright galaxies at lower $z \sim 10–14$ (see, e.g., P. A. Oesch et al. 2016; P. Arrabal Haro et al. 2023; S. L. Finkelstein et al. 2024; S. Carniani et al. 2024b).

Surprisingly, none of the current extragalactic JWST surveys have managed to securely identify brighter objects beyond $z > 16$. No bright ($M_{UV} \lesssim -18$) galaxy candidates have been identified by substantially larger areas covered by CEERS and PRIMER (M. B. Bagley et al. 2023; C. T. Donnan et al. 2023; S. L. Finkelstein et al. 2024). The same is the case for moderately large and deep observations from JADES (D. J. Eisenstein et al. 2023) and NGDEEP (G. C. K. Leung et al. 2023; M. B. Bagley et al. 2024). What our observations seem to indicate is that this dearth of bright galaxies at high $z$ is driven by the shape of the UVLF itself, which can be accommodated by very efficient star formation on top of the steep halo mass function, requiring significantly fainter intrinsic observations to discover $z > 16$ galaxies in reasonable volumes. It is unsurprising, therefore, that the unique sensitivity of GLIMPSE—achieved through a combination of lensing magnification and ultradeep imaging—has managed to potentially reveal some of the earliest galaxies in the Universe.

Finally, by evaluating various star formation scenarios, we demonstrate that $M_{UV} \sim -17$ galaxies at $z \sim 16–17$ can potentially grow into the observed UV-bright galaxies seen by JWST at $z \sim 10–14$, without necessarily violating the standard cosmological framework. The extreme SFE and burstiness coupled with a steep halo mass function are both sufficient to explain the observed trends. This connection further emphasizes the idea that faint, yet numerous, high-$z$ galaxies play a crucial role in the stellar mass assembly of galaxies in the early Universe.

## Acknowledgments

We thank the anonymous referee for their thorough feedback, which helped to improve this manuscript. We thank Steven Finkelstein and Volker Bromm for helpful discussions that helped improve the quality of this manuscript. V.K. acknowledges support from the University of Texas at Austin Cosmic Frontier Center. 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, cofunded by CEA and CNES. The BGU lensing group 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. J.B.M. acknowledges support through NSF grants AST-2307354 and AST-2408637. A.S.-L. acknowledges support from the Knut and Alice Wallenberg Foundation. A.A. acknowledges support by the Swedish
research council Vetenskapsradet (2021-05559). This work has
received funding from the Swiss State Secretariat for Educa-
tion, Research and Innovation (SERI) under contract No.
MB22.00072, as well as from the Swiss National Science
Foundation (SNSF) through project grant 200020_207349.
This work is based on observations made with the NASA/
ESA/CSA James Webb Space Telescope. All of the data
presented in this Letter were obtained from the Mikulski
Archive for Space Telescopes (MAST) at the Space Telescope
Science Institute. The specific observations analyzed can be
accessed via doi: 10.17909/4byn-fe55. These observations are
associated with program #3293. Some of the data products
presented herein were retrieved from the Dawn JWST Archive
(DJA). DJA is an initiative of the Cosmic Dawn Center, which
is funded by the Danish National Research Foundation under

grant DNRF140. 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.

Facilities: JWST, HST.

Software: BAGPIPES (A. C. Carnall et al. 2019), BEAGLE (J. Chevallard & S. Charlot 2016), EAZY (G. B. Brammer et al. 2008), GALFIT (C. Y. Peng et al. 2002), grizli (G. Brammer 2023), sep (K. Barbary 2016), SExtractor (E. Bertin & S. Arnouts 1996).

**Appendix**
**Source Photometry**

In this appendix we present Table 2, which contains the AB
magnitude for each of high-z candidates presented in Figure 3.

<table><caption>**Table 2**
Total AB Magnitude within *D* = 0.2 Circular Aperture</caption><thead><tr><th>ID</th><th>F090W</th><th>F115W</th><th>F150W</th><th>F200W</th><th>F277W</th><th>F356W</th><th>F410M</th><th>F444W</th><th>F480M</th></tr></thead><tbody><tr><td>70467</td><td>&gt;31.90</td><td>&gt;31.90</td><td>&gt;31.70</td><td>&gt;31.80</td><td>30.61 ± 0.19</td><td>30.92 ± 0.25</td><td>&gt;31.10</td><td>30.89 ± 0.28</td><td>&gt;30.10</td></tr><tr><td>72839</td><td>&gt;31.60</td><td>&gt;31.80</td><td>&gt;31.70</td><td>&gt;31.70</td><td>30.55 ± 0.20</td><td>30.86 ± 0.26</td><td>30.53 ± 0.36</td><td>30.87 ± 0.27</td><td>&gt;30.20</td></tr></tbody></table>
ORCID iDs

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

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

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

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

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

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

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

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

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

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

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

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

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

Jeremy Blaizot https://orcid.org/0000-0003-1609-7911

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

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

Gourav Khullar https://orcid.org/0000-0002-3475-7648

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

Ilias Goovaerts https://orcid.org/0009-0007-8470-5946

Michelle Jecmen https://orcid.org/0009-0004-4725-8559

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

Floriane Leclercq https://orcid.org/0000-0002-6085-5073

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

Charlotte Mason https://orcid.org/0000-0002-3407-1785

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

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

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

Erica Nelson https://orcid.org/0000-0002-7524-374X

Joki Rosdahl https://orcid.org/0000-0002-7534-8314

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

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

Maxime Trebitsch https://orcid.org/0000-0002-6849-5375

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

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

References

Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2023, MNRAS, 518, 4755

Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357

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

Atek, H., Chemerynska, I., Wang, B., et al. 2023, MNRAS, 524, 5486

Atek, H., Richard, J., Kneib, J.-P., & Schaerer, D. 2018, MNRAS, 479, 5184

Austin, D., Adams, N., Conselice, C. J., et al. 2023, ApJL, 952, L7

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

Bagley, M. B., Pirzkal, N., Finkelstein, S. L., et al. 2024, ApJL, 965, L6

Balestra, I., Vanzella, E., Rosati, P., et al. 2013, A&A, 559, L9

Barbary, K. 2016, JOSS, 1, 58

Beouchesne, B., Clément, B., Hibon, P., et al. 2024, MNRAS, 527, 3246

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

Bergamini, P., Rosati, P., Mercurio, A., et al. 2019, A&A, 631, A130

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

Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734

Boylan-Kolchin, M. 2023, NatAs, 7, 731

Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils: v1.0.0, Zenodo, doi:10.5281/zenodo.4044744

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

Brammer, G. 2023, grizli, v1.8.2, Zenodo, doi: 10.5281/zenodo.7712834

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

Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

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

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

Calabrò, A., Pentericci, L., Santini, P., et al. 2024, A&A, 690, A290

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

Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2019, MNRAS, 490, 417

Carniani, S., D'Eugenio, F., Ji, X., et al. 2024a, arXiv:2409.20533

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

Casey, C. M., Akins, H. B., Shuntov, M., et al. 2024, ApJ, 965, 98

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

Chabrier, G. 2003, PASP, 115, 763

Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464

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

Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068

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

Conselice, C. J., Adams, N., Harvey, T., et al. 2024, arXiv:2407.14973

Costantin, L., Pérez-González, P. G., Vega-Ferrero, J., et al. 2023, ApJ, 946, 71

Cullen, F., McLure, R. J., McLeod, D. J., et al. 2023, MNRAS, 520, 14

Curtis-Lake, G., Carniani, S., Cameron, A., et al. 2023, arXiv:2306.02465

Dayal, P., Ferrara, A., Dunlop, J. S., & Pacucci, F. 2014, MNRAS, 445, 2545

Dayal, P., Rossi, E. M., Shiralilou, B., et al. 2019, MNRAS, 486, 2336

DeCoursey, C., Egami, E., Pierel, J. D. R., et al. 2025, ApJ, 979, 250

Dekel, A., & Krumholz, M. R. 2013, MNRAS, 432, 455

Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, 523, 3201

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

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

Elíasdóttir, Á., Limousin, M., Richard, J., et al. 2007, arXiv:0710.5636

Endsley, R., Stark, D. P., Lyu, J., et al. 2023, MNRAS, 520, 4609

Endsley, R., Stark, D. P., Whitler, L., et al. 2024, MNRAS, 533, 1111

Feldmann, R., Boylan-Kolchin, M., Bullock, J. S., et al. 2025, MNRAS, 536, 988

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

Ferrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, 522, 3986

Ferrarese, L., Côté, P., Jordán, A., et al. 2006, ApJS, 164, 334

Finkelstein, S. L., Leung, G. C. K., Bagley, M. B., et al. 2024, ApJL, 969, L2
Furtak, L. J., Zitrin, A., Weaver, J. R., et al. 2023, MNRAS, 523, 4568

Gehrels, N. 1986, ApJ, 303, 336

Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279

Gutkin, J., Charlot, S., & Bruzual, G. 2016, MNRAS, 462, 1757

Hainline, K. N., Johnson, B. D., Robertson, B., et al. 2024, ApJ, 964, 71

Harikane, Y., Inoue, A. K., Ellis, R. S., et al. 2025, ApJ, 980, 138

Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, arXiv:2303.11946

Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805

Jaacks, J., Finkelstein, S. L., & Bromm, V. 2018, MNRAS, 475, 3883

Jeffreys, H. 1961, Theory of Probability (3rd edn.; Oxford: Oxford Univ. Press)

Kassiola, A., & Kovner, I. 1993, ApJ, 417, 450

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Kokorev, V., Brammer, G., Fujimoto, S., et al. 2022, ApJS, 263, 38

Leung, G. C. K., Bagley, M. B., Finkelstein, S. L., et al. 2023, ApJL, 954, L46

Li, Z., Dekel, A., Sarkar, K. C., et al. 2024, A&A, 690, A108

Loeb, A., & Furlanetto, S. R. 2013, The First Galaxies in the Universe (Princeton, NJ: Princeton Univ. Press)

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

Maiolino, R., Scholtz, J., Witstok, J., et al. 2024, Natur, 627, 59

Mason, C. A., Trenti, M., & Treu, T. 2015, ApJ, 813, 21

Mauerhofer, V., & Dayal, P. 2023, MNRAS, 526, 2196

McLeod, D. J., Donnan, C. T., McLure, R. J., et al. 2024, MNRAS, 527, 5004

Meena, A. K., Zitrin, A., Jiménez-Teja, Y., et al. 2023, ApJL, 944, L6

Mirocha, J., La Plante, P., & Liu, A. 2021, MNRAS, 507, 3872

Monna, A., Seitz, S., Greisel, N., et al. 2014, MNRAS, 438, 1417

Morishita, T., & Stiavelli, M. 2023, ApJL, 946, L35

Muñoz, J. B., Mirocha, J., Furlanetto, S., & Sabti, N. 2023, MNRAS, 526, L47

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

Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129

Oke, J. B. 1974, ApJS, 27, 21

Pacucci, F., Dayal, P., Harikane, Y., Inoue, A. K., & Loeb, A. 2022, MNRAS, 514, L6

Pascale, M., Frye, B. L., Diego, J., et al. 2022, ApJL, 938, L6

Pei, Y. C. 1992, ApJ, 395, 130

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097
Pérez-González, P. G., Costantin, L., Langeroodi, D., et al. 2023, ApJL, 951, L1

Polletta, M., Tajer, M., Maraschi, L., et al. 2007, ApJ, 663, 81

Ren, K., Trenti, M., & Mutch, S. J. 2018, ApJ, 856, 81

Richard, J., Claeysens, A., Lagattuta, D., et al. 2021, A&A, 646, A83

Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023, PASP, 135, 028001

Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001

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

Rosdahl, J., Blaizot, J., Katz, H., et al. 2022, MNRAS, 515, 2386

Schmidt, M. 1968, ApJ, 151, 393

Schouws, S., Bouwens, R. J., Ormerod, K., et al. 2024, arXiv:2409.20549

Schwarz, G. 1978, AnSta, 6, 461

Sérsic, J. L. 1963, BAAA, 6, 41

Shipley, H. V., Lange-Vagle, D., Marchesini, D., et al. 2018, ApJS, 235, 14

Steidel, C. C., Giavalisco, M., Pettini, M., Dickinson, M., & Adelberger, K. L. 1996, ApJL, 462, L17

Steinhardt, C. L., Jauzac, M., Acebron, A., et al. 2020, ApJS, 247, 64

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

Sun, F., Helton, J. M., Egami, E., et al. 2024, ApJ, 961, 69

Tacchella, S., Diemer, B., Hernquist, L., et al. 2019, MNRAS, 487, 5416

Topping, M. W., Stark, D. P., Senchyna, P., et al. 2024, MNRAS, 529, 3301

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

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

Wilkins, S. M., Lovell, C. C., Irodotou, D., et al. 2024, MNRAS, 527, 7965

Zavala, J. A., Buat, V., Casey, C. M., et al. 2023, ApJL, 943, L9

Zitrin, A., Fabris, A., Merten, J., et al. 2015, ApJ, 801, 44