© The Authors 2026

1. Introduction

The James Webb Space Telescope (JWST; Gardner et al. 2006, 2023) has revolutionized our ability to detect faint galaxies at high redshifts, thanks to its unprecedented sensitivity and wavelength coverage extending beyond 2 μm. In particular, the identification of near-infrared (NIR) dropout sources – objects that drop out of detection in the bluest available NIRCam filters – has emerged as a powerful approach for uncovering previously hidden galaxy populations. These dropouts can trace a wide range of astrophysical phenomena across cosmic time.

The NIR dropouts can signal the presence of high-redshift (z > 10) galaxies, where the dropout is caused by the Lyman break redshifted into the NIR bands. These galaxies offer a unique window into the epoch of reionization and the early formation of stars, black holes, and large-scale structure (e.g., Barkana & Loeb 2001; Naidu et al. 2022; Carniani et al. 2024). Moreover, pre-launch simulations had predicted JWST’s capability to detect galaxies even at ultrahigh redshifts (UHRs; z > 15, e.g., Cowley et al. 2017; Behroozi et al. 2020; Yung et al. 2020), and several candidates have been photometrically identified (e.g., Austin et al. 2023; Harikane et al. 2023; Conselice et al. 2025; Pérez-González et al. 2025; Castellano et al. 2025) and had their brightness compared with theoretical models and predictions (e.g., Ferrara et al. 2023; Matteri et al. 2025a,b). However, no spectroscopic confirmations have yet been achieved (Robertson et al. 2024; Pérez-González et al. 2024b), due to the extreme faintness of these objects and degeneracies in photometric redshift estimates.

At intermediate redshifts (1 ≲ z ≲ 6), NIR dropouts can unveil heavily dust-obscured galaxies, whereby the lack of detection in bluer bands arises from strong attenuation by cosmic dust. These include massive, ultraviolet (UV) and/or optically dark galaxies (e.g., Wang et al. 2016; Enia et al. 2022) as well as the recently identified high-extinction, low-mass (HELM) population (Bisigello et al. 2023). The latter are especially intriguing, as their combination of low stellar mass and high dust content challenges standard dust production models that associate significant dust enrichment with high stellar mass and strong star formation (e.g., Brinchmann et al. 2004; Dayal et al. 2013; Bisigello et al. 2024a). Moreover, other sources such as galaxies with strong line emission at intermediate redshift (Arrabal Haro et al. 2023; Donnan et al. 2023) can reproduce the photometric properties expected for z > 15 galaxies. A prominent example is CEERS-93316, a z = 4.9 dusty star-forming galaxy whose emission-line spectrum mimicked the Lyman break of a z ∼ 16 source (Harikane et al. 2023; Pérez-González et al. 2023). Finally, even cold, substellar Milky Way objects such as brown dwarfs (Holwerda et al. 2024) characterized by strong molecular absorption features can appear as NIR dropouts, mimicking the photometric colors of both dusty galaxies and UHR systems. The overlapping photometric signatures of these diverse populations lead to substantial degeneracies in spectral energy distribution (SED) fitting and complicate the interpretation of NIR dropouts. To disentangle these overlapping populations, robust photometric analysis pipelines must account for a wide range of physical scenarios, including variation in dust attenuation, star formation histories (SFHs), and possible contamination from stars or active galactic nuclei (AGNs). Identifying the most promising candidates for future follow-up requires not only broad and deep multiwavelength photometry, but also careful modeling that evaluates and compares alternative redshift and physical parameter solutions.

In this study, we conduct a targeted search for faint, previously unpublished F200W-dropout sources in the CEERS survey (Finkelstein et al. 2023; Bagley et al. 2023; Yang et al. 2023; Finkelstein et al. 2025). Our goal is to identify and characterize candidates that may belong to one of the extreme populations described above: either dusty galaxies at low to intermediate redshift, or UHR systems at z > 15. To this end, we apply a multi-pronged SED fitting analysis incorporating diverse SFHs and dust laws, and we introduce CosMix, a novel tool designed to stack the photometric signal of similar sources to improve the robustness of the inferred physical parameters. Our framework explicitly accounts for alternative solutions, including intermediate-z strong line emitters and potential AGNs.

The paper is organized as follows. Section 2 provides an overview of the CEERS survey, describing the data exploited in this work. Sect. 3 presents the catalogs used in the analysis and provides relevant information on the extraction of the photometry. Sect. 4 details the selection technique used to obtain our sample of objects and a discussion of their photometric colors. Sect. 5 outlines the SED-fitting procedure adopted on our objects. Sect. 6 delivers the core results of our analysis showing the results of the adopted SED-fitting procedure, as well as computing the UV luminosity function at UHRs. Sect. 7 contains a discussion of our results and Sect. 8 draws the relevant conclusions, giving prospects for future spectroscopic studies that could target our sources.

Throughout this work, we adopt the cosmological parameters yielded by the latest Planck collaboration release (Planck Collaboration VI 2020) and a Kroupa (2001) initial mass function (IMF). We express all magnitudes in the AB system (Oke & Gunn 1983).

2. CEERS data

This work uses imaging data obtained by JWST’s NIRCam as part of the CEERS survey (P.I. S. Finkelstein; Finkelstein et al. 2017, 2023; Bagley et al. 2023; Yang et al. 2023; Finkelstein et al. 2025). CEERS covers ∼90 arcmin² of the Extended Groth Strip (EGS; Davis et al. 2007) field with JWST imaging and spectroscopy through a 77.2 hours Director’s Discretionary Early Release Science Program. CEERS NIRCam observations are available in three short-wavelength bands (F115W, F150W, F200W) and four long-wavelength bands (F277W, F356W, F410M, F444W). Supplementary coverage in the F090W band is provided by the Cycle 1 GO program “The JWST-legacy narrow-band survey of H-alpha and [OIII] emitters in the epoch of reionization” (P.I. E.o Bañados). This program also includes observations in the F470N band, which will not be utilized in this work due to their low sensitivity. CEERS NIRCam pointings were previously targeted by Hubble Space Telescope (HST) observations both with the Advanced Camera for Surveys (ACS; covering F435W, F606W, and F814W) and Wide Field Camera 3 (WFC3; covering F105W, F125W, F140W, and F160W) as part of the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011).

This work utilizes the latest release of CEERS NIRCam imaging data (v.1.0), processed with the JWST Calibration Pipeline (v.1.13.4; Bushouse et al. 2024) and a customized CEERS pipeline for enhanced 1/f noise subtraction and artifact removal, as detailed in Bagley et al. (2023). Compared to earlier releases (v.0.51 or v.0.2), Data Release (DR) v.1.0 provides significant improvements in removing “snowballs” (cosmic-ray-induced visual artifacts, cleaned using James Davies’ snowblind software¹) and “wisps” (stray-light artifacts; see Tacchella et al. 2023; Rieke et al. 2023), along with enhanced handling of saturated and defective pixels. These refinements are crucial, as some artifacts can mimic compact, red galaxies, potentially overlapping with our dropout selection criterion, though visual inspection often suffices for exclusion. NIRCam images were pixel-aligned to the HST data (Koekemoer et al. 2011; Rest et al. 2023), and all HST images were converted to a pixel scale of 0.03”. The Full Width at Half Maximum (FWHM), Point Spread Function (PSF) and 5σ depth of the relevant CEERS NIRCam and HST observations are reported in Table 1. With respect to previous versions, CEERS DR v.1.0 offers improved astrometric alignment across multiband images and superior background subtraction.

CEERS JWST MIRI data are limited (Yang et al. 2023), with the available filters varying based on the pointing considered. In some cases, only two deep short-wavelength filters (F560W and F770W) are included, while in others, six broadband filters spanning from 5.6 μm to 21 μm are available but with shallower depths. However, as we discuss in Sect. 4, our sources are not covered by CEERS MIRI data (see also Fig. 1).

Fig. 1. RGB mosaic of CEERS DR v.1.0 NIRCam data showing the position of our sample of dropouts, which is described in Sect. 4. Sources belonging to the UHR galaxy candidate sample are represented by colored star markers (see Sect. 7.1), while the strong line emitter CURION is depicted with a diamond marker (see Sect. 7.4). All other F200W dropouts are instead indicated by gray circles. The available CEERS DR 0.6 MIRI pointings overlapping with the EGS NIRCam-covered area are shown as red contours.

3. Photometry

To ensure a comprehensive selection of F200W-dropout sources, we rely on two catalogs constructed from CEERS survey data. The first consists of an updated version of the v0.51 CEERS catalog, originally described in Finkelstein et al. (2024, Steven Finkelstein, private communication). The second is the public CEERS ASTRODEEP-JWST catalog² (Merlin et al. 2024). Using both catalogs is essential for our study, as each employs a different photometric extraction technique, meaning some faint high-redshift sources may be missed in one catalog but captured in the other (as discussed in Sect. 4). Both catalogs are optimized for detecting objects (even very faint ones) in the long-wavelength NIRCam bands, a critical feature given our focus on identifying faint, very red objects. Appendix A provides an overview of the key steps involved in creating both catalogs, as well as the methods used for estimating fluxes and photometric uncertainties.

4. Sample selection

We aim to select F200W-dropout galaxies in both the updated CEERS and CEERS ASTRODEEP-JWST catalogs. This selection is specifically designed to include HELM galaxies as well as z > 10 galaxies without restricting the sample to only blue, relatively dust-free objects (see the discussion in Sect. 2.3 of B23), allowing for the selection of UHR candidates with breaks falling between F200W and F277W (or even at longer wavelengths). We hence require our objects to have:

• S/N ≥ 3in NIRCam/F444W;

• S/N ≤ 3 in every NIRCam/HST filter equal or below 2 μm;

• S/N < 3 in the coadded F090W+F115W+F150W+F200W image.

We enforced a signal-to-noise ratio (S/N) of < 2 for all HST and short NIRCam wavelength bands to ensure our selected objects are all F200W dropouts, enforcing a non-detection in the coadded NIRCam short wavelength images to further ensure the dropout nature of our objects. We created the coadded NIRCam short wavelength by a weighted sum of the F090W, F115W, F150W and F200W mosaics pixel by pixel. The coadded image photometry was subsequently measured using the photutils Python package (Bradley et al. 2024), employing circular apertures sized to match the average fiducial aperture for each source used in the updated CEERS catalog’s and CEERS ASTRODEEP’s JWST photometry. We applied a slightly more lenient threshold on the coadded images (S/N < 3), as the summing process includes some of the background noise. The S/N  ≥  3 requirement in F444W corresponds to a detection threshold of AB  ≲  29.6, while the S/N  ≤  2 non-detection criterion applied to all NIRCam/HST filters with λ ≤ 2 μm implies approximate AB magnitude thresholds of F090W  ≳  30.8, F115W  ≳  30.9, F150W  ≳  30.9, F200W  ≳  31.1, F435W  ≳  30.5, F606W  ≳  30.5, F814W  ≳  30.4, F105W  ≳  29.5, F125W  ≳  29.4, F140W  ≳  29.4, and F160W  ≳  29.4.

While our F200W-dropout selection shares similarities with the HELM galaxy selection of B23, a key difference lies in the underlying adopted catalogs and images in the analysis. The HELM galaxies of B23 are identified in the v.0.51 CEERS Photometric Catalog (Finkelstein et al. 2023), which is based on a co-added F277W+F356W detection image. Our work, on the other hand, utilizes the updated CEERS and CEERS ASTRODEEP-JWST catalogs, which include sources extracted using the coadded F356W+F444W image as a detection image. This approach enables us to identify F444W-detected F200W-dropout sources that were not captured in previous CEERS catalogs. To focus on these novel objects, we exclude any sources appearing in the latest, publicly available release of the CEERS catalog (v0.51.4) within a 0.3″search radius. This procedure grants that there are not overlap between our sample of sources and the updated CEERS HELM catalog from Bisigello et al. (2025), obtained via private communication, which is based on a co-added F277W+F356W detection image covering the total CEERS area.

Our selection criteria leave us with 636 (in the updated CEERS catalog) and 941 (in the CEERS ASTRODEEP-JWST catalog) sources. After a careful visual inspection to remove diffraction spikes, residual artifacts, noise peaks in the F444W band and objects with positions falling in the chip gap of all NIRCam short-wavelength observations, our final sample of F200W dropouts includes 8 updated CEERS objects and 3 CEERS ASTRODEEP-JWST objects, totaling 11 objects. 5 of these 11 F200W-dropout objects lack coverage in either the F090W or F115W bands. We list the IDs, coordinates and names of our F200W-dropout sample in Table 2 and their photometry in Table C.1.

4.1. Cross-match with other catalogs

We further checked that our sample of sources has no overlap with the CEERS little red dot (LRD) sample found in Kocevski et al. (2025) and Taylor et al. (2025), or the high-z sample presented in Finkelstein et al. (2024)³, or the CEERS UHR sample presented in Yan et al. (2023, highlighted via dropout selections; see their Sect. 5.1), or the sample of HST-dark galaxies revealed by Barrufet et al. (2023) in CEERS. We checked for potential spectroscopic counterparts of the sources in v0.5 of the CEERS NIRCam F356W Wide Field Spectroscopy (WFSS) catalog covering pointings 5, 7, 8 and 9 and CEERS DR 0.7 NIRSpec catalog, finding that our objects are not covered by spectroscopic observations (no objects were discarded during this step). Moreover, despite the few objects falling near the edges of the available CEERS MIRI pointings, no MIRI observations are available for our sources (with the only exception of “CURION”, i.e., U-112842; as we discuss in Sect. 7.4). A few of our sources are also covered by (sub)millimeter data. Source A-76468 falls near the low-sensitivity edge of a proprietary NOEMA 1.1 mm map (NOEMA program W20CK, PIs V. Buat and J. Zavala, obtained via private communication), and it is undetected at a 3σ level (with a 3σ upper limit of 0.9 mJy). Sources VICIADGO (U-34120), LIZZAN (U-75985), NASPE (U-80918) and CURION (U-112842) are instead covered by 850 μm SCUBA-2 observations (described in Zavala et al. 2017 and Geach et al. 2017), and all of these are undetected as well (with 3σ upper limits of ∼3.4 mJy/beam). We checked during our SED-fitting procedure that the inclusion of these (sub)millimeter upper limits in the available photometry does not change our results in any way. To determine the physical properties of these objects we therefore rely on broadband JWST/NIRCam, HST/ACS, and HST/WFC3 photometry, with the aim of corroborating the present analysis with deep MIRI imaging, wider and deeper NOEMA observations and/or future spectroscopic data targeting our sample.

It is important to highlight that, since both the updated CEERS catalog and the CEERS ASTRODEEP-JWST catalog target the same field, a large overlap between them is expected, with the two catalogs having 63 034 objects in common (exploiting a search radius of 0.3″). As for our sample, we find that all three CEERS ASTRODEEP-JWST sources are also present in the updated CEERS catalog. This means that the updated CEERS catalog’s photometry of these sources did not meet our selection criteria due to the different techniques exploited to compute photometry errors in the two catalogs (see Sect. 3), with the updated CEERS catalog’s ones being higher on average. An example of this is illustrated in Fig. 2, which compares the updated CEERS catalog and CEERS ASTRODEEP-JWST photometric errors of common objects between the two catalogs in a reference band (F356W). The aim of this is to show how the updated CEERS catalog’s photometric errors are systematically greater than CEERS ASTRODEEP-JWST ones, especially for higher error values. We do not delve into pinpointing the details behind these discrepancies in photometric uncertainty estimations, as it falls outside the scope of our paper.

Fig. 2. Updated CEERS catalog photometry errors vs. CEERS ASTRODEEP-JWST photometry errors for the sources in common between the two catalogs (displayed as gray dots) for a sample band (F356W). The bisector of the plot is represented as a dashed black line surrounded by a 1σ confidence interval (shaded black area) representing the typical dispersion of the data points around the bisector. A power-law fit is displayed to showcase how the updated CEERS catalog’s errors tend to be higher with respect to CEERS ASTRODEEP-JWST ones.

Additionally, we confirmed that not all updated CEERS catalog sources in our sample are found in the CEERS ASTRODEEP-JWST catalog. This is the case for U-53105, which is absent from the CEERS ASTRODEEP-JWST catalog within a 0.5″search radius from its updated CEERS catalog’s position. This emphasizes the value of combining multiple catalogs with different source extraction methods to capture a broader array of sources.

4.2. Comparison to red color-selected objects

This section summarizes the key color-magnitude properties of our F200W-dropout sample, subsumed in Figs. 3, 4, and 5. Our sources are rather heterogeneously distributed in the [F277W – F444W] diagram (Fig. 3), with five of them lying within the 95% CEERS source density contours, and only three of them with [F277W – F444W] > 2.5. Four sources (ARCERIO, NASPE, VICIADGO and LIZZAN) fall within the average errors below the 5 σ depth limit in the F444W band, making their selection particularly uncertain. However, we still attempt to characterize their physical properties to enhance our preliminary results with future, deeper observations targeting this field. Notably, four of our objects (CURION, U-53105, A-76468 and U-100588) are within 1σ of the color-magnitude space utilized in Pérez-González et al. (2024a) to select LRDs in the SMILES/MIRI survey (i.e., [F277W – F444W] > 1 and F444W < 28, plus the additional requirement of F150W–F200W< 0.5). Additionally, a fifth, object (U-106373) is at a 2σ-distance from the LRDs color-magnitude selection space. We assess the possibility of these objects being LRDs in Sect. 5 during our SED fitting runs.

Fig. 3. [F277W – F444W] color-magnitude diagram for our F200W dropouts. The F200W-dropout sample’s objects are shown as colored stars (UHR candidates), pea green diamond (CURION, a strong line emitter), cyan cross (stacked UHR galaxies) or gray circles (all other objects in the sample). For CURION, U-100588 and U-106373 we report 1σ upper limits. Black dots in the background represent the colors of CEERS galaxies, with shaded gray contours enclosing 50% (inner), 80% (middle), and 95% (outer) of the total CEERS source density. The F444W 5σ depth is marked as a dashed black line, while the shaded brown area represents the color-magnitude space occupied by LRDs in Pérez-González et al. (2024a, albeit such a selection also requires F150W–F200W < 0.5). Gray error bars represent the average errors for objects with colors that are not upper limits.

Fig. 4. [F150W – F444W] color-magnitude diagram for our F200W dropouts. The F200W-dropout representation scheme is the same adopted in Fig. 3, as well as the F444W 5σ depth threshold and x axis average errors. In addition, the HST-dark galaxy selection technique from Gottumukkala et al. (2024) is represented by a continuous brown line.

Fig. 5. [F277W–F356W] versus [F200W–F277W] color-color diagram for our sources. The F200W-dropout representation scheme is the same adopted in Figs. 3 and 4, as well as the x axis average errors. The black rectangle highlights the UHR LBG selection for 15 < z < 20 galaxies by Castellano et al. (2025), whereas the shaded blue area corresponds to the 15 < z < 20 LBG selection adopted in Kokorev et al. (2025).

The [F150W – F444W] color-magnitude diagram (Fig. 4) expresses rather well the extreme properties of our sources. In this color-magnitude space, all our sources lay outside the contour encompassing the 95% of CEERS sources density adopting 1σ upper limits, suggesting that they are outliers in terms of [F150W – F444W] colors with respect to the bulk of other CEERS sources (note that our objects are all undetected in the F150W band due to the adopted selection technique). Moreover, approximately half of our sources (CURION, BUDIARA, U-100588, U-106373, A-76468) exhibit colors falling above the selection threshold for HST-dark galaxies adopted in (i.e., F150W - F444W > 2.1; Gottumukkala et al. 2024), further emphasizing the remarkable color properties of our objects.

Finally, a [F277W – F356W] versus [F200W – F277W] color-color diagram (Fig. 5) proves useful for evaluating which of our objects exhibit colors consistent with those expected for 15 < z < 20 Lyman-break galaxies. Typically, 9 < z < 15 Lyman-break galaxies are identified as red [F115W – F444W] objects with flat SED continua redward of the break (see, e.g., Castellano et al. 2022). For comparing our sources to the expected colors of 15 < z < 20 Lyman-break galaxies, we employ the selection adopted in Castellano et al. (2025) to select such objects. This color-magnitude selection space is illustrated in Fig. 5 as a black rectangle. We also reported as a shaded blue area the color-magnitude space corresponding to the selection adopted in Kokorev et al. (2025) to select 15 < z < 20 galaxies. This area corresponds to requesting [F200W – F277W] > 1.0, [F200W – F277W] > 1.2 + 2.0[F277W – F356W], [F277W – F356W] < 0.5. Remarkably, five of our objects (VICIADGO, ARCERIO, U-53105, U-106373 and A-26130) are consistent within 1σ (i.e., considering the average errors) with the selection of Kokorev et al. (2025) and seven of them with the reddened selection exploited in Castellano et al. (2022), indicating that their colors align with those of UHR galaxies in the 15 < z < 20 range. Moreover, as reported in Fig. 3 of Holwerda et al. (2024), the F200W-F277W color of BDs in CEERS is clustered around 0 – the color distribution of our objects would hence tend to exclude that they are BDs in the Milky Way.

5. SED-fitting procedure

To characterize the physical properties of our objects, we fit their available photometry using different SED-fitting codes, and compared their results. First, we performed SED fits with Bagpipes (Carnall et al. 2018), relying on the Bruzual & Charlot (2003) stellar population synthesis (SPS) model. For each object, we run separate SED-fitting runs exploiting all the available parametric SFHs in Bagpipes – an exponentially declining SFH (Mortlock et al. 2017; Wu et al. 2018; McLure et al. 2018), a delayed exponentially declining SFH (Ciesla et al. 2017; Chevallard et al. 2019), a log-normal SFH (Diemer et al. 2017; Cohn 2018), and a double power-law SFH (Ciesla et al. 2017; Carnall et al. 2018). As is discussed in Carnall et al. (2019), some of these models may be biased against very high-z solutions due to their inherent inability to reproduce rising SFHs (especially the exponential and the delayed exponential models which trace declining SFHs; see e.g. Reddy et al. 2012; Harikane et al. 2022). However, we are interested in investigating the impact of different SFHs on the final results of our SED-fitting runs, and we discuss this matter further in Sect. 6. Note, however, that the stellar mass depends marginally on the shape of the parametric form of the adopted SFH, and we expect the mass estimates between different SFHs to be broadly consistent (see, e.g., Santini et al. 2015; Ciesla et al. 2024). We also performed our Bagpipes runs testing two different dust attenuation laws: a Small Magellanic Cloud (SMC) dust attenuation law (Prevot et al. 1984; Bouchet et al. 1985) and a standard Calzetti law (Calzetti et al. 2000) to compare the results. To enhance sensitivity to strong line emitter solutions, which may occupy a narrow region in the probability distribution of galaxies’ physical parameters, we increase the number of Bagpipes’s live points (i.e., the walkers used in the program’s Markov Chain Monte Carlo sampling exploited by MultiNest algorithm; Buchner et al. 2014; Feroz et al. 2019) from the default 400 to 2000 (A. Carnall, private communication), ensuring a robust exploration of this parameter space. Additionally, we include a nebular emission component in all Bagpipes model fits (modeled through the CLOUDY photoionization code; Ferland et al. 2017), allowing the ionization parameter to reach logU = −1. A detailed list of the model fit parameters and their prior ranges is provided in Table B.1. To account for systematic errors, we add 5% of the flux estimates in quadrature to the photometric uncertainties for each run.

Since Bagpipes does not natively account for an AGN component, we also fit our objects with CIGALE (Boquien et al. 2019), which instead allows one to model an AGN contribution to the galaxies’ overall emission. Since the aim of these CIGALE runs is to reveal potential lower-z AGN contaminants, we did not adopt a wide set of SFHs as was done for our Bagpipes runs, and we limited ourselves to only exploring the effect of a double exponential SFH. Furthermore, our CIGALE SED-fitting runs rely on the single stellar populations library of Bruzual & Charlot (2003), including a nebular emission component with CLOUDY and a modified Charlot & Fall (2000) dust attenuation law. Finally, we include in the CIGALE model fit template AGN emission models by Fritz et al. (2006), allowing AGN fractions, which correspond to the AGN-heated dust contribution to the total infrared luminosity (integrated over 1–1000 μm), to vary between 0% and 75%. A complete list of the free parameters of our CIGALE fit and the related value grids are available in Table B.2. For each run, we accounted for systematic errors by adding 5% of the flux estimates in quadrature to the photometric uncertainties.

Finally, we carried out SED fitting using the Python implementation of EAZY⁴ (Brammer et al. 2008), employing the following broad set of SED-fitting templates:

• The EAZY template set from Brammer et al. (2008), based on the input set of Grazian et al. (2006) following the Blanton & Roweis (2007) algorithm plus a template accounting for dusty galaxies;

• Default fsps templates based on flexible SPS models (specifically, we adopted the tweak_fsps _QSF_12 _v3.param template);

• sfhz template set, which incorporates redshift-dependent SFHs, excluding those that begin earlier than the Universe’s age at any given epoch, along with a realistic emission-line model at z ∼ 8 by Carnall et al. (2023, specifically, we tested for every object the agn_blue_sfhz_13.param, blue_sfhz_13.param, carnall_sfhz_13.param and corr_sfhz_13.param templates);

• High-redshift optimized templates⁵ (hot; Steinhardt et al. 2023), a custom set of two synthetic templates optimized to reproduce the expected astrophysics of galaxies at high redshifts. These templates assume hotter gas temperatures of 45K for redshifts z = 8 − 12 and 60K for z > 12 in star-forming regions, resulting in bottom-lighter or top-heavier IMFs.

Our EAZY fits were conducted between z = 0.01 and z = 25 with a step of 0.01.

6. Results

6.1. SED-fitting results

Our Bagpipes analysis highlights five sources in our F200W-dropout sample with a significant P(z) volume at UHRs (z > 15; see Table 3): BUDIARA (U-31863), VICIADGO (U-34120), LIZZAN (U-75985), NASPE (U-80918) and ARCERIO (A-22691). For every SFH and dust extinction model assumed, these objects always showcase bimodal P(z)s, featuring one HELM-like solution at z < 8 (i.e., a dusty, low-mass galaxy solution) and a UHR solution with negligible dust, meaning that the available JWST/NIRCam, HST/ACS and HST/WFC3 photometry alone cannot describe exact results of their nature. Due to the bimodal nature of their P(z)s, we did not simply adopt the best-fit values output by Bagpipes as the definitive properties of our objects. This is because these values would reflect the physical characteristics of a median solution between the low- and high-redshift peaks of the bimodal P(z) distributions. Instead, we divided each object’s P(z) into two separate components: one for z ≤ 8 and another for z > 8. For each of these redshift regimes, we determined the best-fit solutions maximizing their likelihood. This method allows us to rigorously account for both possible interpretations of each object’s nature. Consequently, we report two distinct sets of best-fit physical parameters for each source: one representing the low-redshift, z≤ 8 solution and the other the high-redshift, z > 8 solution. The likelihood maximization for the z ≤ 8 and the z > 8 solutions draws from the different assumed SFHs, which means that the solutions at z ≤ 8 and z > 8 for the same object could be associated with two different SFHs. We report our fiducial Bagpipes best fit SEDs and normalized z ≤ 8 and z > 8 P(z)s obtained following this procedure in Fig. D.1, with the inferred physical parameters for both solutions being listed in Appendix D.

Another highlight from our analysis is CURION (U-112842). For this particular object we were able to include in our analysis data from the Cycle 2 Proposal MEGA Mass Assembly at Cosmic Noon: MIRI EGS Galaxy and AGN Survey Kirkpatrick et al. (2023); P.I. A. Kirkpatrick, obtained via private communication). These MIRI observations target the EGS field, covering ERS NIRCam observations in the F770W, F1000W, F1500W, and F2100W bands. The best-fit SED produced by our Bagpipes runs including MIRI data favors a log-normal SFH, with a best-fit redshift of z ∼ 5.4, an intermediate dust content (A_v ∼ 0.87), a low mass (logM/M_⊙ ∼ 7) and extreme line emission, with logU ∼ −3.9, approaching the lower limits of the priors assigned to this parameter in the fits. All the best-fit parameters inferred by these SED-fitting runs are listed in Appendix F. As shown in Fig. 7, CURION is consistent within 1σ with the main sequence of star-forming galaxies. Fits yielded by EAZY and CIGALE do not hint at the presence of an AGN, favoring dusty solutions at z < 10, yielding, however, best-fit chi-squared values up to 7 or 8 times the ones returned by Bagpipes in the z ∼ 5.5 strong line emitter scenario. Finally, since CURION is unresolved in both the F356W and F444W bands (see Appendix G), we explored the possibility of it being a substellar Milky Way object such as a BD. To do so, we fit its photometry with L and T dwarf models from Burrows et al. (2006), spanning different metallicities (−0.5 ≤ [Fe/H] ≤ 0.5), effective temperatures (700 K ≤T ≤ 2300 K) and surface gravity accelerations (10^4.5cm s⁻² ≤ g ≤ 10^5.5cm s⁻²). Our BD template fits return a best-fit chi-squared value of χ_BD² ∼ 20.96, up to 6 times higher than the Bagpipes one.

Bagpipes fits for the source U-53105 yield a primary low-redshift peak at z ∼ 2.6 (adopting an exponential SFH, which yields the best chi-squared value for this object), with a low-probability P(z) up to z ∼ 15, while a compatible main redshift solution is obtained by varying the adopted SFH. A log-normal SFH maintains the primary peak at z < 3, while raising the P(z) volume predicted at z > 5, with two strong secondary peaks at z ∼ 9 and z ∼ 15. However, EAZY indicates a z < 5 solution for this object, even if some templates assign a strong probability beyond z > 20, with the blue_sfhz, carnall _sfhz, corr _sfhz, hot _45k and hot _60k all yielding ∫_z > 15dz P(z)≥0.95. A similar trend is followed by CIGALE’s P(z), showcasing a marked bimodal distribution, with principal P(z) peaks falling at z < 8 and z > 15, not indicating the presence of an AGN. However, we note that the high solution volume at z > 20 predicted by EAZY and CIGALE is not reproduced by our Bagpipes runs even when adopting a log-normal SFH. Furthermore, U-53105 appears resolved in all NIRCam long-wavelength bands (see Appendix G), ruling out the possibility of it being a BD.

Source A-26130’s Bagpipes fits adopting delayed, exponential, and double power-law SFHs suggest a best-fit logM/M_⊙ ∼ 8.78, with fairly high dust-obscuration (A_v ∼ 3), located at z ∼ 2. While these fits also reveal a high-redshift tail, the solution volume for such scenarios remains minimal. In contrast, a log-normal SFH fit implies a slightly higher mass logM/M_⊙ ∼ 8.7, retaining a primary peak at z = 2 with the addition of a secondary 15 < z < 20 peak. The EAZY P(z) further supports the low-redshift solution but also assigns a primary peak at z ∼ 20 for the majority of the adopted templates, backed up by CIGALE, which does not reveal a significant contribution by AGN emission. The source is resolved in the F277W band (see Appendix G), and we can exclude it from being a BD.

The analysis of source A-76468 with Bagpipes indicates a primary peak at z = 3, with a stellar mass of logM/M ⊙ ∼9 and a dust attenuation of Av ∼ 4. These results are consistent across all tested SFHs. When a log-normal SFH is adopted, the fits also reveal a more extended high-redshift tail, along with secondary peaks at z = 8 and z = 15, reproduced by our CIGALE runs. The EAZYP(z) complements these findings by displaying a narrow primary peak at z ∼ 7 and a secondary peak at z = 3. This source is resolved in both the F277W and F356W bands, excluding the possibility of it being a BD. We report the best-fit Bagpipes SEDs of U-53105, A-26130, and A-76468 in Fig. D.3.

Unfortunately, source U-106373 yields a rather unconstrained P(z)s, and we refrain from reporting its best-fit physical parameters. Bagpipes fits of its photometry, using delayed, exponential, and double power-law SFHs, suggest that this is a relatively massive (9 < logM/M_⊙ < 10) and dusty (A_v ∼ 4) source at redshift z∼5, with a secondary 10 < z < 15 P(z) peak. Adopting a log-normal SFH introduces a tertiary P(z) peak at z > 20. However, the A_v posteriors for all SFHs, while showing a weak probability peak, remain generally unconstrained. Additionally, both EAZY and CIGALE yield rather scattered P(z)s – more data are needed to further constrain the nature of this source.

Finally, source U-100588 is detected in only two NIRCam bands, making it a F356W dropout. This object occupies the edge of the UHR 15 < z < 20 candidate color-magnitude selection space (see Fig. 5), and was characterized in a dedicated work (Gandolfi et al. 2026) with the help of a ∼0.8h-exposure NIRSpec spectrum made available by the CAPERS collaboration.

6.2. UV luminosity function at z ∼ 17

In this section we compute the UV luminosity function of our UHR candidates under the assumption that their true redshift corresponds to the high-redshift solution (z > 8). Given that the best-fit photo-z values in this scenario are mutually consistent within their uncertainties, we assume that all sources lie within the same redshift bin. We then divide our sample into two absolute UV magnitude bins: −20.5 < M_UV < −19.5 and −19.5 < M_UV < −18.5. Given that LIZZAN has a M_UV ∼ −18.491 (i.e., falling outside the faintest bin) as well as having the lowest P(z) volume at z > 8 (see Table 3), we remove it from our computations, and therefore use the remaining four sources. We then compute the luminosity function using the standard 1/V_max method (Avni & Bahcall 1980), adopting the statistical formalism described in Gehrels (1986, see our Table 5). The resulting luminosity function is shown in Fig. 6, alongside points from other works that have attempted to extend the UVLF measurements to similar redshifts (Castellano et al. 2025; Pérez-González et al. 2025; Weibel et al. 2025; Whitler et al. 2025). Overall, we find good agreement with values obtained by comparable studies in CEERS and other deep extragalactic fields.

Fig. 6. UV luminosity function estimate from our work (blue stars), alongside estimates from Castellano et al. (2025, green diamonds), Pérez-González et al. (2025, yellow circles), Weibel et al. (2025, black triangles, representing upper limits), and Whitler et al. (2025, red triangles, again representing upper limits). The best-fit power law is represented as a dashed black line.

Due to the narrow range in absolute UV magnitudes probed by our sample, and the relatively large associated uncertainties, a double power-law fit is not well constrained. For this reason, we adopt a single power-law model to describe the data (i.e., log ϕ = αM_UV + β, a trending similar to the one followed by the faint end of the UV luminosity function). However, this approach comes with a caveat: caution should be exercised when integrating the resulting luminosity function over a much broader magnitude range than directly observed. In fact, the single-slope approximation may not reliably capture the true shape of the luminosity function at both the bright and faint ends. That said, the observed data are well fit by a single power law characterized by a slope α ∼ 0.4, roughly corresponding to a double power-law faint-end slope as α = −(α_DBL/0.4+1), meaning that α_DBL ∼ −2.18 (i.e., roughly the same faint-end slope retrieved in Donnan et al. 2024 at z > 9). The best-fit power-law intercept is β ∼ 4.39. Assuming M_* = −19.88 (i.e., the M_* value reported by Donnan et al. 2024 at z = 9.7), we estimate the equivalent normalization ϕ_* of a DPL as log(ϕ_*) = αM_* + β = −4.95 Mpc⁻³ mag⁻¹, between the estimates normalizations by Donnan et al. 2024 (−4.74 Mpc⁻³ mag⁻¹) and Weibel et al. 2025 (< −5.05 Mpc⁻³ mag⁻¹) at z = 17.

7. Discussion

In this section we discuss our results, weighing the possible solutions for each candidate examined through our SED fitting procedure.

7.1. Five candidates at z> 15

The z > 8 solution for all of our five UHR candidates consistently favors a log-normal SFH. This is not surprising, since log-normal SFHs may reproduce both rising and declining star formations, while the other models may be less flexible in representing the SFH shapes found in the early Universe. All the z > 8 best-fits redshifts of our sources are included in the 17 < z < 18 range, and are accompanied by negligible dust, with the highest A_v best-fit value being the one of ARCERIO (A_v ∼ 0.27), while in general other candidates are less dusty (A_v < 0.2). The predicted best-fit metallicities are rather unconstrained, showing errors > 60% of the measurements, indicating that further data are needed in order to obtain reliable estimates in this sense. The best-fit values for the ionization parameter logU are compatible with the expected SED’s flat continuum at these redshifts, with no strong line emission. We note that the best-fit galaxy age parameters defining the SFH shape are quite unconstrained and degenerate with each other – this is a consequence of the limited number of available detections in our dataset, and we preferred not to assume any particular recipe to constrain these values and bias our results. However, our results clearly favor SFHs that started within the first Gyr of the Universe’s life. Finally, we checked the consistency of the z > 8 best-fit stellar mass of each object with Λ cold dark matter (ΛCDM) in terms of baryon-to-star conversion efficiency in Sect. 7.3.

The z ≤ 8 solutions of our sample of UHR candidates tend to favor HELM-like solutions, similar to the ones obtained for objects discussed in Sect. 7.5. The inferred typical best-fit masses amount to logM/M_⊙ < 8, high dust content (A_v > 4) and z < 2. The main difference between our sample and the HELMs of B23 is that the majority of our objects in the z≤ 8 solution are dusty enough to all be within 1.5σ from the upper limit we set for the dust attenuation index (A_v = 6) – our objects could therefore be a particularly dusty class of HELMs, which is a rather intriguing possibility. We thus performed a comparison between the z≤ 8 solution of our UHR candidates and the CEERS HELM galaxy sample of (Bisigello et al. 2025); encompassing 4141 sources) in terms of stellar masses versus dust attenuation (see Fig. 10). The A_v measurements’ centroids for all our UHR candidates in their z ≤ 8 solutions fall on the edge or beyond the contour enclosing 95% of the CEERS HELM galaxies. In the z ≤ 8 scenario, our objects could be probing an extreme dust extinction regime if compared to other HELM galaxies given their stellar mass, which could enhance our understanding of how cosmic dust is produced in such low-mass systems. This regime is similar to the one probed by the z ∼ 4.0 solution of the NIRCam-dark “Cerberus” galaxy (Pérez-González et al. 2024b). If the z ≤ 8 nature of our UHR candidates was to be confirmed, their intrinsic faintness would not be due to their extreme distance, but it would rather be attributable to their exceptional dust content. However, we stress that the lack of data in MIR and (sub)millimeter wavelengths or spectroscopic observations makes the dust content characterization of these sources particularly uncertain. We shall therefore conduct further studies of their dust properties as soon as other data become available.

The only exception to this is constituted by LIZZAN, whose z ≤ 8 solution is located at a higher redshift (z ∼ 4.31, albeit with a rather large lower error bound) and has a smaller dust attenuation (A_v ∼ 2.04) with respect to the other UHR candidates. LIZZAN’s best-fit z ≤ 8 solution’s SED indeed shows prominent emission lines, making LIZZAN a potential strong line emitter contaminant, and suggesting that our set-up is capable of correctly accounting for strong line emitter solutions for our galaxy sample. The middle and bottom panels of Fig. 7 highlight that in the z < 8 solutions our UHR candidates are all consistent within 1σ with the main sequence of star-forming galaxies.

Fig. 7. Comparison between galaxies in our sample and the main sequence of star-forming galaxies. The upper panel displays four of our UHR candidates assuming their best-fit z ≤ 8 solutions (represented as colored stars) alongside their median stack (represented as a cyan cross). Since the best-fit z < 8 solution for LIZZAN is at a higher redshift (z ∼ 4.31) with respect to its companion UHR candidates, we checked its consistency with the main sequence in a separate plot (bottom panel) alongside CURION (portrayed as a pea green diamond), falling in a similar redshift window. The main sequence is modeled after the parametrization of Popesso et al. (2023) and is represented in each panel as a black line computed at the average redshift between the displayed objects (⟨z⟩∼1.00 and ⟨z⟩∼4.89, respectively), alongside a 1σ confidence interval (shaded gray area).

We repeated our Bagpipes SED-fitting runs adopting a Calzetti dust attenuation law to assess the impact of different extinction assumptions on the inferred properties for our UHR candidate sample. We found that altering the dust attenuation law does not significantly change our results, yielding consistent physical parameters for our objects. The agreement between the two dust attenuation laws suggests that, in both the z ≤ 8 and z > 8 scenarios, our UHR candidates lack a significant 2100 Å bump. Such a bump is linked to the presence of carbon lattice grains in the interstellar medium (ISM), and is the key difference between the SMC and Calzetti dust models (Hensley & Draine 2023). Indeed, some parameters inferred from our SED-fit analysis are intrinsically poorly constrained due to the limited number of photometric data points, leading to large uncertainties. Both the Calzetti and SMC attenuation laws yield similarly bimodal solutions for the UHR candidates studied here, though they differ slightly in how the solution space in the P(z) plane is distributed between HELM and UHR probabilities. Notably, the Calzetti law appears to favor UHR solutions in certain objects, even when using SFHs such as delayed or exponentially delayed models, further emphasizing the intriguing nature of our UHR candidates. Since the physical properties derived using the Calzetti law align with those obtained using the SMC law, we adopt the SMC law as the fiducial attenuation model for the remainder of our analysis. We plan to revisit the impact of different dust attenuation laws on the inferred properties of our galaxies as additional data at longer wavelengths become available, which should better highlight subtle distinctions between the two.

Our EAZY runs exploiting the wide set of templates listed in Sect. 5 tend to favor, on average, the z > 8 solutions of our UHR galaxy candidates’ samples in terms of the inferred P(z)s. Both BUDIARA and VICIADGO feature an integral probability ∫_z > 15dz P(z) > 0.99 for all the templates used in the analysis, indicating a strong preference for z > 15 best-fit solutions. LIZZAN shows an integral probability ∫_z > 15dz P(z) > 0.99 for the majority of the EAZY templates utilized in our runs, with the exception of hot-45k templates (∫_z > 15dz P(z) > 0.68) and hot-60k (∫_z > 15dz P(z) > 0.79). These latter two templates seem to allocate part of LIZZAN’s P(z) volume to a secondary solution at z < 10, somewhat compatible with the secondary, lower-z strong line emitter solution predicted by Bagpipes for this source. However, both the hot-45k and hot-60k still predict that the bulk of the P(z) volume for LIZZAN is located at z > 15, favoring the UHR nature of this source. EAZY templates predict a ∫_z > 15dz P(z) > 0.99 for NASPE, with the exception of the fsps templates, yielding ∫_z > 15dz P(z)∼10⁻³. The P(z) predicted for NASPE by this particular template reproduces the one yielded by Bagpipes at z < 8, backing the HELM galaxy scenario for this source and assigning a negligible solution volume at z > 8. This makes NASPE the only source in our five UHR galaxy sample to have a strongly bimodal EAZY P(z). Finally, our EAZY runs predict for ARCERIO 0.78 < ∫_z > 15dz P(z) < 0.95, depending on the template used. This source thus shows a mild P(z) bimodality, with the bulk of the EAZY P(z)s always lying beyond z = 15, regardless of the template used.

Our CIGALE runs systematically back up the ultrahigh-z nature of our sources, with P(z) peaks aligning to the ones predicting by other SED-fitting codes (with the exception of ARCERIO, predicting an even higher best-fit redshift of ∼21.2). Furthermore, our CIGALE runs did not reveal any significant AGN component associated with our UHR sample, yielding best-fit predictions for both the AGN fraction always ≤0.4 and, in many cases, compatible with zero.

Finally, we remark that all our UHR candidates are resolved in at least one NIRCam band (see Appendix G). This is enough to exclude the possibility of them possibly being substellar Milky Way objects such as BDs.

7.2. Stacking the UHR galaxy candidate sample with CosMix

Despite the inherent degeneracy introduced by the bimodality of their P(z) distributions, the objects in our UHR galaxy sample exhibit remarkably consistent physical properties, allowing us to treat them as a homogeneous group for further analysis. To improve the S/N of our five UHR candidates and to mitigate the degeneracy between the z ≤ 8 and z > 8 solutions, we chose to stack the sample. This approach allows us to analyze their combined median photometric properties, providing a more robust framework for understanding their nature.

To stack the sources in our UHR candidate sample, we used CosMix (Gandolfi et al., in prep.), a new Python-based stacking tool designed to combine the multiwavelength signal of faint sources with homogeneous physical properties and enhance their characterization. CosMix will soon be available as a free web application, running on a dedicated machine hosting a multiwavelength set of data spanning different fields available for stacking. CosMix extracts image cutouts from input catalogs (ID, RA, Dec), combines them, and outputs downloadable FITS files of the desired size. It also computes stacked photometry with customizable settings, producing an SED and a FITS file with errors and upper limits, ready for fitting with any preferred program. We report the salient passages of our setup for stacking in Appendix E. After obtaining the median stacked photometry through CosMix, we applied the same SED-fitting procedures described in Sect. 5 for fitting single objects. The median UHR candidates stack shows a S/N < 3 in all bands blueward of F277W. Some residuals are visible in F200W due to the narrow contrast scale of the cutout, with a S/N ∼ 1.7, not enough to be associated with a detection. Furthermore, the F115W stack has a S/N ∼ 2.8 – this value, which is not high enough to be associated with a sound detection, is likely attributable to noise accumulation during the stacking process, which could not be mitigated due to the limited number of objects in our sample. As visible in Fig. 5, the color-magnitude properties of the stacked UHR candidates are broadly in line with the ones of the single objects in the sample, being consistent with both the UHR galaxies’ selections reported in the plot.

Our Bagpipes stacked photometry fits highlight that the stacking procedure was not able to completely remove the P(z) bimodality of our UHR candidates, which is expected due to the narrow wavelength range of the available data. However, we now find that even delayed and delayed exponential SFHs admit a non-negligible solution volume at z > 8 (as visible by looking at the gray P(z) in Fig. 8), marking an important difference with respect to single-object fits. A strong improvement is also visible by looking at the median P(z) yielded by our EAZY fits, which now shows a completely negligible solution volume at z < 10, with ∫_z > 15dz P(z)∼0.99 for all the templates used in the analysis. Our CIGALE fits yield consistent best-fit physical properties with respect to the other codes. These results further corroborate the intriguing nature of our UHR candidates and their potential as promising targets for future spectroscopic studies. To infer the physical properties of our stacked object, we split its Bagpipes P(z)s in two as we did with our single sources. Finally, the physical parameters inferred in this way (listed in Table E.2) are consistent with the z≤ 8 and z > 8 solutions of our single UHR candidates. As we did with single objects, we assessed the stack’s best-fit stellar mass compatibility with ΛCDM in terms of baryon-to-star conversion efficiency in Sect. 7.3. We find that the stacked UHR galaxy sample is consistent with the main sequence of star-forming galaxies in the z ≤ 8 solution (see the top panel of Fig. 7).

Fig. 8. Best-fit SEDs, P(z)s, and NIRCam cutouts for the median stacked UHR candidate sample, obtained exploiting CosMix. The plot’s color scheme is the same followed in Fig. D.1.

7.3. Consistency with ΛCDM

To corroborate whether the z > 8 solutions of our five UHR candidates could be realistic within the standard ΛCDM framework, we checked the consistency of their best-fit stellar masses and baryon-to-stars conversion efficiencies with the ΛCDM model. The underlying idea behind this test is simple: the hierarchical formation of cosmic structures predicted by ΛCDM sets a limit on the maximum dark matter halo mass expected at a given redshift. The mass of dark matter halos (M_halo) can be in turn related to the stellar mass (M_*) of their host galaxies as M_* = ϵf_bM_halo, where ϵ is the efficiency of baryon conversion into stars, and f_b = 0.158 (Planck Collaboration VI 2020) is the cosmological baryon fraction. Assuming a certain dark matter halo mass function, one can thus calculate a theoretical stellar mass limit at each redshift.

In Fig. 9 we show the best-fit z > 8 UHR candidates’ (and their median stack’s) stellar mass estimates versus their best-fit redshift in the z > 8 solution. The black lines in the plot represent theoretical stellar mass limits at each redshift assuming different baryon-to-star conversion efficiencies ϵ, whereas the gray area in the plot represents unrealistic values for the baryon-to-star conversion efficiency ϵ > 1 – a galaxy falling in this area would indicate a tension with ΛCDM. We computed the different ϵ value curves assuming a Sheth & Tormen (1999) dark matter halo mass function and integrating it to determine the mass above which one halo is expected within CEERS’ survey volume, assuming sliding redshift bins of Δz = 2 (i.e., following the approach outlined in Leung et al. 2025). During this process, we made use of the Python COLOSSUS⁶ package (Diemer 2018) to deal with cosmological calculations.

Fig. 9. Best-fit stellar masses of our UHR candidates (and their median stack) in their z > 8 solutions versus best-fit redshift. Each candidate is represented as a colored star, with the median stack being depicted as a cyan cross. The dash-dotted, dotted, dashed, and continuous black lines represent different values of the baryon-to-star conversion efficiency (0.14, 0.2, 0.75, and 1, respectively), while the shaded gray area corresponds to nonphysical efficiency values > 100%.

Our analysis shows that all the UHR candidates in their z > 8 solutions are well below the ϵ = 1 threshold within 1σ, indicating that they are consistent with the ΛCDM framework. BUDIARA and ARCERIO require slightly higher efficiencies with respect to other objects in the sample in order to match the inferred stellar mass, being consistent within 1σ with the ϵ = 0.2 curve, implying higher efficiency values with respect to the ones inferred for z < 4 galaxies (ϵ ∼ 0.14, see Chworowsky et al. 2024). Other well-characterized high-redshift galaxies such as PENNAR (Rodighiero et al. 2023; z_photom = 12.1), GHZ2 (Castellano et al. 2024; z_spec = 12.34), or JADES-GS-z14-0 (Carniani et al. 2024, z_spec = 14.32) yield ϵ values around ϵ ∼ 0.3. These high star-formation efficiencies would also be in line with recent COSMOS-Web results showing an increased value for ϵ in z > 5 galaxies (Shuntov et al. 2025). Other UHR candidates’ z > 8 solutions (plus the stack) showcase instead a milder conversion of baryons into stars with ϵ ≲ 0.1, similar to the efficiency attained at a lower redshift by JADES-GS-z14-1 (Carniani et al. 2024, z_spec = 13.90). All in all, this test indicates that the high-redshift solutions of our objects are plausible within the ΛCDM framework.

As a further check, we estimated the physical sizes of our UHR candidates in all long-wavelength NIRCam bands assuming their z > 8 solutions (see Appendix G). We then compared the inferred sizes with the empirical effective radius versus redshift relation found in Westcott et al. (2025), extrapolated beyond z = 17 (see Fig. G.1). All our UHR candidates in their z > 8 turned out to be consistent with the extrapolated Westcott et al. (2025) relation within 2σ (with the exception of LIZZAN in the F277W and F356W bands, which instead falls within 2σ from the aforementioned relation).

7.4. “CURION” – A strong line emitter at z∼ 5

With a best-fit Bagpipes redshift of z ∼ 5.5, CURION is consistent with being a strong line emitter galaxy contaminant, underscoring the important of multiwavelength data in UHR galaxy searches. In fact, preliminary EAZY SED-fitting runs relying on JWST/NIRCam, HST/ACS and HST/WFC3-only photometry yielded a primary peak in CURION’s P(z) at z ∼ 13, with a secondary solution at z ∼ 9. Bagpipes runs were able to recover these solutions assuming a log-normal SFH, while for more conservative SFHs the primary P(z) peak turned out to be at z ∼ 5.5. However, the z ∼ 5.5 solutions predicted by Bagpipes relying on JWST/NIRCam, HST/ACS and HST/WFC3 photometry alone always yielded extremely high dust attenuation index values, compatible with the adopted parameter’s prior upper limit (A_v ≥ 6). Even if CURION appears undetected in the available MEGA MIRI data, the inclusion of photometric information between 7.7 and 21 μm in the SED-fitting procedure drastically changes this picture, helping to constrain CURION’s dust content and determine its redshift more robustly. Finally BD fits yield to the conclusion that CURION is indeed a galaxy and likely not a substellar Milky Way object.

7.5. Three dusty dwarf galaxies

The best-fit physical properties derived for U-53105 make us conclude that U-53105 is in all likelihood a dusty galaxy at z ∼ 2.6. A-26130 indeed falls within the 15 < z < 20 UHR selection space (see Fig. 5), yet its P(z) does not exhibit the distinctly bimodal distribution typically associated with other UHR candidates discussed in this work. This lack of a pronounced high-redshift solution might stem from the absence of a significant spectral break in the best-fit inferred SED, consistent with the faint emission observed in the F200W band for this source. The faint visibility of the source in the F200W band also suggests that the high-redshift solution found by EAZY and CIGALE is likely spurious. Preliminary evidence therefore suggests that A-26130 is likely a dusty galaxy at z ∼ 2 rather than a bona fide UHR source in the 15 < z < 20 range. Finally, for A-76468 both Bagippes, EAZY and CIGALE do not support the interpretation of this object as a UHR galaxy – the high stellar masses (logM/M_⊙ > 9.5) predicted by the high-redshift solutions lead us to favor the conclusion that this source is a dusty galaxy located at z = 3. All three sources exhibit physical properties similar to the ones of other HELM galaxies in CEERS, tracing their high-stellar mass end.

We show that U-53105, A-26130 and A-76468 are all compatible within 1σ with the main sequence of star-forming galaxies in Fig. 11. Furthermore, we detailed a comparison between the HELM galaxies discussed in this section with the CEERS HELM galaxy sample of Bisigello et al. (2025) in terms of stellar masses versus dust attenuation (see Fig. 10). Our three HELM galaxy candidates fall at more than 3σ from the contour enclosing 50% of the CEERS HELM galaxies in the stellar mass versus dust attenuation plane, indicating that they belong to the more massive end of the HELM galaxies distribution. In particular, A-26130 and A-76468 show higher best-fit stellar masses with respect to U-53105, and they lie within 1σ from the 95% contour. Investigating this massive end of the HELM galaxies distribution could be crucial to characterize how stellar masses and dust attenuation correlate for these objects, highlighting potential differences with the expected mass versus dust attenuation relation for the other galaxies in the CEERS sample. However, data at longer wavelengths (such as in the MIR or (sub)millimeter bands) or spectroscopic follow-ups would be needed to reduce the uncertainties affecting estimates of the dust content of the three HELM galaxies presented here.

Fig. 10. Comparison between the 4141 sources in the CEERS HELM galaxy sample by Bisigello et al. (2025) and the HELM solutions retrieved in this work in the form of a logarithmic stellar mass versus dust attenuation index plot. Blue contours represent isodensity lines in the distribution of the updated CEERS HELM galaxy sample by Bisigello et al. (2025). The dashed, dash-dotted, and dotted black lines are enclosing, respectively, 50%, 80%, and 95% of the total HELM distribution. Colored star markers represent our UHR candidates adopting their z ≤ 8 solutions, with their median stack being represented by a cyan cross. The three colored circles instead represent our strongest HELM candidates discussed in Sect. 7.5. Finally, the brown hexagon represents the best-fit stellar mass and dust attenuation index for the NIRCam-dark “Cerberus” galaxy (Pérez-González et al. 2024b) in its z ∼ 4 solution (fitting with the synthesizer-AGN code and adopting a double-exponential SFH with a Calzetti dust attenuation law).

Fig. 11. Comparison between U-53105, A-26130 and A-76468 to the main sequence of star-forming galaxies, represented as purple circles. The main sequence is modeled after the parametrization of Popesso et al. (2023) and is represented in each panel as a black line computed at the average redshift between the displayed objects (⟨z⟩∼2.76 alongside a 1σ confidence interval (shaded gray area).

8. Conclusions

Our analysis revealed a sample of faint, previously missed F200W dropouts in the CEERS survey. These objects are interesting candidates for future observations aimed either at soundly characterizing their SED in the MIR/(sub)millimeter regime via the addition of new photometric information (thus extending their presently available wavelength coverage beyond the F444W band), or at spectroscopic follow-up. We investigated the physical properties of these sources, mostly relying on the available NIRCam photometric data and complementing them with MIR/(sub)millimeter data when available. Here, we recap the most significant results of our analysis:

• We revealed three z < 4 HELM galaxies (U-53105, U-26130, and A-76468; i.e., dwarfs characterized by a puzzlingly high dust content) that occupy the high-mass end of the stellar mass distribution of HELM objects in CEERS. This regime can crucially help in characterizing how the stellar mass of these galaxies correlates with their dust content.

• We highlighted one (CURION, U-112842) strong line emitter at z ∼ 5. Preliminary EAZY runs exploiting NIRCam-only photometry yielded a primary z ∼ 13 solution, but the inclusion of MIRI upper limits from the MEGA collaboration revealed its true nature.

• We reported a sample of five faint objects with bimodal P(z)s – our analysis suggests that these are either extremely dusty HELM galaxies or UHR objects between 15 < z < 20, whose stellar masses are not in tension with ΛCDM cosmology. Either possibility is extremely intriguing – yet, in order to be further thoroughly investigated, more longer-wavelength photometric data or detailed spectroscopic follow-ups are needed.

• We outlined an analysis pipeline that combines different SED-fitting codes, templates, SFHs, dust extinction models, and stacking that could be helpful in future UHR object searches, helping to discriminate between known lower-z contaminants and assess the P(z) solution volume of each source.

We plan to extend our set of multiwavelength photometric data describing the SEDs of the objects in the sample analyzed in this work. In particular, our object CURION highlights the discriminating power of MIR data in the SED-fitting procedure of these dropout objects, even when these data are just upper limits. Moreover, deeper observations in the far-infrared/(sub)millimeter bands could help in constraining the dust content of our sources, both from current facilities such as NOEMA or future ones such as NASA’s PRobe far-Infrared Mission for Astrophysics (PRIMA; see e.g., Bisigello et al. 2024b). Extending the multiwavelength coverage of these sources would be pivotal to test the impact of different and more detailed dust extinction models on their best-fit solutions, nonparametric SFHs, and/or different SPS models. On the other hand, we also plan to refine our analysis including new SED-fitting codes such as GalaPy (Ronconi et al. 2024) implementing a two-component dust absorption and/or re-radiation model that is unbiased against the physics of dust grains and integrated with high-resolution stellar populations for the study of primordial galaxies.

Spectroscopic coverage of our sources would also be pivotal to soundly constrain their nature. A promising opportunity in this sense is constituted by the potential inclusion of our sources in the upcoming CAPERS program. A few of our objects (the five UHR candidates) are likely too faint to probe their continuum emission, yet it would still be possible to observe or rule out emission lines from solutions at z < 10, casting rather strong constraints on their true nature.

To date, only a very small number of HELM galaxies have been spectroscopically confirmed (Bisigello et al. 2024a; Castellano et al. 2025). A larger sample of spectroscopically confirmed HELM galaxies is important to better characterize dust production mechanisms in such low-mass galaxies, especially if, as is hinted at in our analysis, some of the targets have such an extreme dust content given their stellar mass with respect to other similar objects. Moreover, our analysis suggests that HELM galaxies can mimic the color-color properties of UHR, 15 < z < 20 galaxies (see Fig. 5) – an in-depth knowledge of their spectral properties will allow us to filter them out effectively from future UHR object searches, building a comprehensive set of templates.

On the other hand, a few of our objects (the ones discussed in Sect. 7.1) show a significant solution volume beyond z > 15. Although our analysis is not sufficient to draw robust and definitive conclusions on the UHR nature of these sources, future follow-ups could enhance the UHR solution volumes for these galaxies, contributing to breaking the z = 15 galaxy barrier. Characterizing the statistics of galaxies beyond z = 15 will have deep astrophysical implications, contributing to understanding how the first galaxies and their black holes formed and evolved. However, the discovery and characterization of UHR galaxies also has deep astroparticle implications. As has been shown by preliminary works using early release JWST data, the characterization of UHR sources falling near the rest-frame UV luminosity function’s knee may help in constraining the particle properties of dark matter, to the point of successfully discriminating between competing paradigms beyond the CDM framework (Gandolfi et al. 2022).

Finally, the pipeline developed in this work to characterize the physical properties of faint, F200W-dropout objects in CEERS will be applied in other fields covered by deeper JWST observations. Finding other candidates with bimodal P(z)s would have the benefit of enlarging our stacked UHR candidate sample to further constrain the properties of these objects, while potentially outlining promising UHR candidates.

Acknowledgments

G. G, G. R. and B. V. are supported by the European Union – NextGenerationEU RFF M4C2 1.1 PRIN 2022 project 2022ZSL4BL INSIGHT. M. G. acknowledges support from INAF under the following funding schemes: Large Grant 2022 (project “MeerKAT and LOFAR Team up: a Unique Radio Window on Galaxy/AGN co-Evolution”) and Large GO 2024 (project “MeerKAT and Euclid Team up: Exploring the galaxy-halo connection at cosmic noon”). We thank Adam Carnall for the stimulating discussions and the suggestions to refine our Bagpipes SED-fitting procedure. We thank D. Giacobazzi for suggesting the names of our UHR galaxy candidate sample. Finally, we thank the anonymous referee for the helpful suggestions and comments, which contributed in improving this manuscript.

References

Appendix A: Catalogs

Here, we describe the catalogs exploited in our analysis: the updated CEERS catalog and the ASTRODEEP-JWST catalog, hilighting key differences.

A.1. The updated CEERS catalog

The photometry process on which this catalog is based broadly follows what is described in Finkelstein et al. (2024), and we summarize here the salient passages. This catalog, which relies on the recent CEERS DR 1.0, was built by PSF-matching ACS and NIRCam images blueward than F277W to the F277W image, while computing appropriate correction factors for images with larger PSFs. Photometry was then performed with Source Extractor (SExtractor; Bertin & Arnouts 1996; v2.25.0) in two image mode, using the co-added F277W+F356W native resolution image as the detection image. Co-addition was performed exploiting a weighted sum, using pipeline-generated error (ERR) map weights (i.e., error maps including Poisson noise) for the measurement image and 1/sqrt(WHT) weights (i.e., effective RMS maps not including Poisson noise) for the measurement one, so that Poisson noise, included in the ERR maps, does not affect detections. A first hot-mode run (i.e., a run optimized for the detection of faint objects; see, e.g., Galametz et al. 2013) is performed setting DETECT_THRESH = 1.2, DETECT_MINAREA = 5, DEBLEND_NTHRESH = 32 and DEBLEND_MINCONT = 0.0001, with these values derived iteratively via visual inspection to minimize spurious sources while maximizing completeness. Then, a second cold-mode run (i.e., optimized to extract and deblend efficiently the brightest sources) is performed adopting DETECT_THRESH = 2.4, DETECT_MINAREA = 50, DEBLEND_NTHRESH = 32 and DEBLEND_MINCONT = 0.01. Both runs exploit a 4-pixel-wide top-hat convolution kernel with a kernel size of 5×5 pixels. Fiducial fluxes are measured in small Kron apertures (Kron factor = 1.1 and minimum radius = 1.6, as in Finkelstein et al. 2023), applying two aperture corrections – one based on the ratio of the default larger MAG_AUTO Kron aperture to our small aperture and the second to compensate for the flux missed on larger scales, derived by source injection simulations. Flux uncertainties were estimated by randomly placing non-overlapping apertures, with diameters ranging from 0.1” to 1.5”, in empty areas of the image, avoiding the segmentation map. The normalized median absolute deviation was then fit exploiting a four-parameter function, repeating the fit for each aperture as a function of their area. The flux error for each object was determined using this fit based on the area of the object’s aperture. To build the final catalog, cold mode objects were added to hot objects not included in the cold catalog. Moreover, a third run exploiting F444W observations as detection image was added to the two pre-existing catalogs. The final updated CEERS catalog consists of 62,752 hot run sources, 30,703 cold run sources and 11,326 F444W sources, for a total of 104,781 objects.

A.2. The CEERS ASTRODEEP-JWST catalog

The other part of our sample of objects relies on the CEERS ASTRODEEP-JWST catalog (Merlin et al. 2024). ASTRODEEP-JWST gathers data from eight JWST NIRCam observational programs (GLASS-JWST, UNCOVER, DDT2756, GO3990, PRIMER, JADES, NGDEEP and CEERS DR 0.5 and 0.6), exploiting carefully chosen detection parameters designed to maximize the detection of high-redshift faint extended sources. To build this catalog, the astrometries of JWST and HST images were re-matched down to the milli-arcsec level. The RMS and ERR maps were then adjusted by applying a correction factor, which was determined by comparing the flux dispersion of a set of simulated point-like sources (generated using WebbPSF) with their nominal errors. Science images were then smoothed with a Gaussian convolution filter with FWHM = 0.14” (close to the F356W and F444W FWHMs) to optimize for lower surface brightness features. Due to this smoothing, the CEERS ASTRODEEP-JWST catalog may be more sensitive to extended sources, while the updated CEERS one is more sensitive to point sources. Photometry was then performed via SExtractor (v2.8.6) adopting a weighted stack of the convolved F356W and F444W images as detection image. The SExtractor run parameters closely resemble the ones used for the CANDELS campaign (Koekemoer et al. 2011; Grogin et al. 2011), adopting DETECT_MINAREA = 5, DEBLEND_NTHRESH = 32 and DEBLEND_MINCONT = 0.0001 (see Table 3 of Merlin et al. 2024 for further details). The final CEERS ASTRODEEP-JWST catalog amounts to 82,547 objects (flagged as detections with a S/N> 2).

Appendix B: SED-fitting priors

We report here the priors and grids utilized during our SED-fitting runs with Bagpipes and CIGALE.

Appendix C: Photometry

We report here the photometry of the F200W dropouts found in our analysis in all the available HST (top table) and NIRCam (bottom table) bands.

Appendix D: Physical properties of the UHR galaxy candidate sample

We report in this section the inferred physical properties yielded by our fiducial Bagpipes runs assuming a SMC dust attenuation law for both the z > 8 (top table) and z≤ 8 (bottom table) solutions, alongside the inferred best-fit χ² for both cases.

Fig. D.1. Best-fit SEDs and redshift probability distributions P(z) obtained with the procedure described in Sect. 5 for the UHR galaxy candidates and their median stack. The top section of each panel shows 1.5″×1.5″ NIRCam cutouts of the source. The left inset displays the observed photometry (black circles) with 1σ upper-limit uncertainties (black triangles). The Bagpipes best-fit SED at z > 8 is shown in red, with its 0.1 dex confidence interval (shaded red area), while the best-fit SED at z ≤ 8 – obtained independently over the same full prior range – is shown in gray. The right inset reports the normalized redshift probability distributions from different SED-fitting codes: the P(z) of the Bagpipes model that minimizes χ2 at z > 8, shown in red over the full prior range (0 ≤ z ≤ 25); the CIGALE posterior (dashed yellow line); and the mean EAZY posterior averaged over all templates described in Sect. 5 (purple). The vertical dashed black line marks the z = 8 threshold.

Fig. D.1. Continued.

Fig. D.1. Continued.

Fig. D.2. Best-fit SEDs, P(z)s and NIRCam cutouts for CURION. The plot’s color scheme is the same followed in Fig. D.1.

Fig. D.3. Best-fit SEDs and P(z)s for U-53105 (top), A-26130 (middle) and A-76468 (bottom). The top section of each plot shows 1.5” × 1.5” NIRCam cutouts for the source. The left inset displays the observed photometry (black circles) with 1σ upper-limit uncertainties (black triangles). The best-fit SED derived with BAGPIPES is shown as a red curve with a 0.1 dex confidence interval (shaded red area). The right inset reports the normalized redshift probability distributions from different SED-fitting codes: the P(z) of Bagpipes in red; the CIGALE posterior in yellow; and the mean EAZY posterior averaged over all templates described in Sect. 5 (purple). The vertical dashed black line marks the z = 8 threshold.

Fig. D.3. Continued.

Appendix E: Stacked cutouts

This appendix lists the salient passages we followed while stacking our UHR candidate sample. We used CosMix to create 2.7” × 2.7” (i.e., 90 × 90 pixels) median cutouts for our UHR candidate sample. We then degraded the PSF of such stacked cutouts to match the F444W’s PSF – this filter, the reddest available, has the largest PSF with a FWHM of 0.161”. All NIRCam filters from F090W to F444W and the three HST/ACS filters (F435W, F606W, and F814W) have smaller FWHMs than F444W, meaning they will be PSF-matched (see Table 1). However, the HST/WFC3 filters have larger FWHMs and thus do not require convolution. PSFs for each filter were obtained by identifying PSF stars in the CEERS images, following the procedure described by Finkelstein et al. (2022). We then used the PYPHER⁷ Python routine (Boucaud et al. 2016) to create convolution kernels and match the PSFs of all filters to F444W. The images were then convolved with their respective kernels to ensure consistency across all bands. After obtaining the PSF-matched stacked cutouts, we needed to select an appropriate photometry aperture size. We determined the aperture size by fitting a Gaussian function to the radial flux profiles for each band featuring a detection (see Fig. E.1). For each band, we used a circular region centered on the stacked source with a radius equal to one third of the cutout size, ensuring that the source was well contained within this area. We then measured the average pixel value outside this circular region to estimate the average stacked background level. Using this background estimate, we interpolated each Gaussian profile to find the radius where the flux matched the background level, repeating the process for each band. We then averaged these aperture sizes across the different bands to obtain a fiducial aperture size passing this value to CosMix to compute the stacked photometry, finding a fiducial value of 0.21” (i.e., 7 pixels). We then used CosMix to perform aperture photometry on the chosen aperture in each stacked band. Stacked photometry errors were estimated using a bootstrap technique, which involves random resampling with replacement from the measured fluxes to create multiple synthetic datasets. The median flux is calculated for each bootstrap sample, resulting in a distribution of medians. The overall median of these medians estimates the central tendency, while the standard deviation provides an approximation of the uncertainty in the measured flux. This method effectively captures variability and accommodates outliers or non-normal distributions in the data, yielding a robust error estimate. We summarize the so-obtained photometry in Table E.1, while SED-fitting results are reported in Table E.2 adopting a SMC dust attenuation law. Note that not all objects are covered in all bands. During stacking, we excluded from SED-fitting any stacked cutouts created from fewer than 2/3 of the total objects (e.g., the F105W band, which has limited coverage in CEERS), though we included the relative flux measurements in the stacked SED plot for completeness.

Fig. E.1. Signal in the convolved long-wavelength NIRCam cutouts is shown for our stacked UHR candidate sample, with insets corresponding to F277W (top left), F356W (top right), F410M (bottom left), and F444W (bottom right). The blue line represents the signal, in μJy, for each pixel in the stacked image as a function of its distance from the center of the cutout (measured in pixels). A dashed green line depicts the Gaussian fit used to estimate the aperture photometry radius, while a dashed red line indicates the mean background level. The vertical dashed blue lines mark the inferred radius for fiducial aperture photometry.

Appendix F: Other F200W dropouts

We list here the physical properties inferred for CURION (U-112842) and the three HELM galaxies (U-53105, A-26130 and A-76468) found in our analysis.

Appendix G: Estimating the sizes of our galaxy candidates

In this section we report estimates of the sizes of our five UHR galaxy candidates (BUIDARA, VICIADGO, LIZZAN, NASPE, ARCERIO), our three HELM candidates (U-53105, A-26130 and A-76468) and CURION. We estimated their effective radii by fitting their two-dimensional flux distribution in each NIRCam long-wavelength broadband channel with a Sérsic profile, given by

$$\begin{array}{c}\hfill I(R)=I_{e}exp\left\{q_s[{\left(\frac{R}{R_e}\right)}^{1/n_s}-1]\right\},\end{array}$$

with I_e being the intensity calculated at the effective radius R_e (i.e., the radius enclosing half of the total light of the source). The overall shape of the profile is determined by the axis ratio q_s and the Sérsic index n_s. To perform the fit, we utilized the code galight⁸ (Ding et al. 2020). This code provides estimations of the angular sizes of our objects in arcseconds in each band, as well as the best-fit values of the Sérsic profile’s axis ratio q_s and index n_s, and estimates of the effective radii in kiloparsecs upon the assumption of a reference redshift. During our galight runs, we set the code’s detection threshold to S/N = 1 in order to fit the shape of our sources despite their intrinsic faintness.

We report the inferred size values in Table G.1 (F277W band), Table G.2 (F356W band) and Table G.3 (F444W band) for our UHR candidates adopting their z > 8 best-fit redshifts. We found that all our UHR candidates are resolved in the F277W band (except for LIZZAN), meaning that these objects are unlikely BDs (which would appear unresolved in all NIRCam bands). The only resolved source in the F356W band is VICIADGO, while the only resolved source in the F444W band is ARCERIO. Our objects seem to have consistent physical sizes (within the measurement errors) across the different bands. Deeper observations would be needed to enhance this comparison, since the intrinsic faintness of our sources translates into large errors in the effective radii estimates. The best-fit Sérsic index values we retrieved are generally < 1, possibly indicating that our UHR candidates are galaxies with peculiar morphologies rather than compact spheroids (see, e.g., Fig. 13 in Ormerod et al. 2023). However, we caveat that our measured Sérsic indices could be also explained as a surface brightness dimming effect on a galaxy with an intrinsic sersic index of ∼ 1 (see, e.g., Westcott et al. 2025).

In the three panels of Fig. G.1, we checked the consistency of the effective radii of our UHR candidates in the F277W, F356W, and F444W bands with the effective radius versus redshift relation found in Westcott et al. (2025). Such a relation was obtained by fitting the effective radii distribution of a sample of 521 galaxies between 6.5 < z < 12.5 from different JWST-covered fields (CEERS, JADES GOOD-S, NGDEEP, SMACS0723, GLASS and PEARLS), and it reads as

$$\begin{array}{c}\hfill R_e=2.74\pm 0.49{(1+z)}^{-0.79\pm 0.08}.\end{array}$$

Fig. G.1. Effective radii (in kiloparsecs) vs. redshift for our UHR candidate sample in the F277W, F356W, and F444W bands, respectively. We report their physical sizes estimated with galight vs. their z > 8 best-fit redshift solutions as colored stars. As a comparison, we report the empirical effective radius vs. redshift scaling retrieved in Westcott et al. (2025) as a continuous black line with 1σ (gray) and 2σ (light gray) confidence intervals.

To enable a comparison with our UHR candidates, we extrapolated this relation to redshifts beyond z ∼ 17 – we stress that this redshift regime is not yet probed by spectroscopically confirmed candidates, and should be thus treated as indicative. We found that all of our UHR galaxy candidates in their z > 8 solutions are consistent within 2σ with the relation by Westcott et al. (2025) extrapolated at UHRs (with the only exception of LIZZAN, which is at a 3σ distance from the Westcott et al. 2025 relation in both the F277W and F356W bands). Our analysis does not highlight any major inconsistency between our inferred best-fit UHR galaxy sizes and the extrapolated relation by Westcott et al. (2025).

We repeated our size estimation with galight for our three HELM candidates and CURION, listing the obtained results in Table G.4, Table G.5 and Table G.6. Our galight run was not able to converge while fitting the size of CURION in the F277W filter due to the extreme faintness of the object in this band; while all other sources are resolved in the F277W band. Instead, only U-53105 and A-76468 are resolved in the F356W band, whereas only U-53105 is resolved in the F444W band. Our angular size estimates seem to back up the hypothesis that these objects are not BDs, albeit we do not have this confirmation for CURION. Similarly to our UHR candidates, these objects showcase compatible physical sizes across the different bands considered.

Finally, in the three panels of Fig. G.2, we perform a comparison between the best-fit sizes for our objects and the relation of Westcott et al. (2025) in an identical fashion to Fig. G.1. All of our HELM candidates are consistent within 2σ with the Westcott et al. (2025) relation, while CURION is consistent with it within 3σ. As in the case of our UHR candidates, our analysis does not highlight any major inconsistency between our inferred best-fit galaxy sizes and the extrapolated relation by Westcott et al. (2025).

Fig. G.2. Same as Fig. G.1, but for our three HELM candidates and CURION.

All Tables

All Figures

Fig. 1. RGB mosaic of CEERS DR v.1.0 NIRCam data showing the position of our sample of dropouts, which is described in Sect. 4. Sources belonging to the UHR galaxy candidate sample are represented by colored star markers (see Sect. 7.1), while the strong line emitter CURION is depicted with a diamond marker (see Sect. 7.4). All other F200W dropouts are instead indicated by gray circles. The available CEERS DR 0.6 MIRI pointings overlapping with the EGS NIRCam-covered area are shown as red contours.

In the text

Fig. 2. Updated CEERS catalog photometry errors vs. CEERS ASTRODEEP-JWST photometry errors for the sources in common between the two catalogs (displayed as gray dots) for a sample band (F356W). The bisector of the plot is represented as a dashed black line surrounded by a 1σ confidence interval (shaded black area) representing the typical dispersion of the data points around the bisector. A power-law fit is displayed to showcase how the updated CEERS catalog’s errors tend to be higher with respect to CEERS ASTRODEEP-JWST ones.

In the text

Fig. 3. [F277W – F444W] color-magnitude diagram for our F200W dropouts. The F200W-dropout sample’s objects are shown as colored stars (UHR candidates), pea green diamond (CURION, a strong line emitter), cyan cross (stacked UHR galaxies) or gray circles (all other objects in the sample). For CURION, U-100588 and U-106373 we report 1σ upper limits. Black dots in the background represent the colors of CEERS galaxies, with shaded gray contours enclosing 50% (inner), 80% (middle), and 95% (outer) of the total CEERS source density. The F444W 5σ depth is marked as a dashed black line, while the shaded brown area represents the color-magnitude space occupied by LRDs in Pérez-González et al. (2024a, albeit such a selection also requires F150W–F200W < 0.5). Gray error bars represent the average errors for objects with colors that are not upper limits.

In the text

	Fig. 4. [F150W – F444W] color-magnitude diagram for our F200W dropouts. The F200W-dropout representation scheme is the same adopted in Fig. 3, as well as the F444W 5σ depth threshold and x axis average errors. In addition, the HST-dark galaxy selection technique from Gottumukkala et al. (2024) is represented by a continuous brown line.
In the text

Fig. 5. [F277W–F356W] versus [F200W–F277W] color-color diagram for our sources. The F200W-dropout representation scheme is the same adopted in Figs. 3 and 4, as well as the x axis average errors. The black rectangle highlights the UHR LBG selection for 15 < z < 20 galaxies by Castellano et al. (2025), whereas the shaded blue area corresponds to the 15 < z < 20 LBG selection adopted in Kokorev et al. (2025).

In the text

	Fig. 6. UV luminosity function estimate from our work (blue stars), alongside estimates from Castellano et al. (2025, green diamonds), Pérez-González et al. (2025, yellow circles), Weibel et al. (2025, black triangles, representing upper limits), and Whitler et al. (2025, red triangles, again representing upper limits). The best-fit power law is represented as a dashed black line.
In the text

Fig. 7. Comparison between galaxies in our sample and the main sequence of star-forming galaxies. The upper panel displays four of our UHR candidates assuming their best-fit z ≤ 8 solutions (represented as colored stars) alongside their median stack (represented as a cyan cross). Since the best-fit z < 8 solution for LIZZAN is at a higher redshift (z ∼ 4.31) with respect to its companion UHR candidates, we checked its consistency with the main sequence in a separate plot (bottom panel) alongside CURION (portrayed as a pea green diamond), falling in a similar redshift window. The main sequence is modeled after the parametrization of Popesso et al. (2023) and is represented in each panel as a black line computed at the average redshift between the displayed objects (⟨z⟩∼1.00 and ⟨z⟩∼4.89, respectively), alongside a 1σ confidence interval (shaded gray area).

In the text

	Fig. 8. Best-fit SEDs, P(z)s, and NIRCam cutouts for the median stacked UHR candidate sample, obtained exploiting CosMix. The plot’s color scheme is the same followed in Fig. D.1.
In the text

Fig. 9. Best-fit stellar masses of our UHR candidates (and their median stack) in their z > 8 solutions versus best-fit redshift. Each candidate is represented as a colored star, with the median stack being depicted as a cyan cross. The dash-dotted, dotted, dashed, and continuous black lines represent different values of the baryon-to-star conversion efficiency (0.14, 0.2, 0.75, and 1, respectively), while the shaded gray area corresponds to nonphysical efficiency values > 100%.

In the text

Fig. 10. Comparison between the 4141 sources in the CEERS HELM galaxy sample by Bisigello et al. (2025) and the HELM solutions retrieved in this work in the form of a logarithmic stellar mass versus dust attenuation index plot. Blue contours represent isodensity lines in the distribution of the updated CEERS HELM galaxy sample by Bisigello et al. (2025). The dashed, dash-dotted, and dotted black lines are enclosing, respectively, 50%, 80%, and 95% of the total HELM distribution. Colored star markers represent our UHR candidates adopting their z ≤ 8 solutions, with their median stack being represented by a cyan cross. The three colored circles instead represent our strongest HELM candidates discussed in Sect. 7.5. Finally, the brown hexagon represents the best-fit stellar mass and dust attenuation index for the NIRCam-dark “Cerberus” galaxy (Pérez-González et al. 2024b) in its z ∼ 4 solution (fitting with the synthesizer-AGN code and adopting a double-exponential SFH with a Calzetti dust attenuation law).

In the text

	Fig. 11. Comparison between U-53105, A-26130 and A-76468 to the main sequence of star-forming galaxies, represented as purple circles. The main sequence is modeled after the parametrization of Popesso et al. (2023) and is represented in each panel as a black line computed at the average redshift between the displayed objects (⟨z⟩∼2.76 alongside a 1σ confidence interval (shaded gray area).
In the text

Fig. D.1. Best-fit SEDs and redshift probability distributions P(z) obtained with the procedure described in Sect. 5 for the UHR galaxy candidates and their median stack. The top section of each panel shows 1.5″×1.5″ NIRCam cutouts of the source. The left inset displays the observed photometry (black circles) with 1σ upper-limit uncertainties (black triangles). The Bagpipes best-fit SED at z > 8 is shown in red, with its 0.1 dex confidence interval (shaded red area), while the best-fit SED at z ≤ 8 – obtained independently over the same full prior range – is shown in gray. The right inset reports the normalized redshift probability distributions from different SED-fitting codes: the P(z) of the Bagpipes model that minimizes χ2 at z > 8, shown in red over the full prior range (0 ≤ z ≤ 25); the CIGALE posterior (dashed yellow line); and the mean EAZY posterior averaged over all templates described in Sect. 5 (purple). The vertical dashed black line marks the z = 8 threshold.

In the text

	Fig. D.1. Continued.
In the text

	Fig. D.1. Continued.
In the text

	Fig. D.2. Best-fit SEDs, P(z)s and NIRCam cutouts for CURION. The plot’s color scheme is the same followed in Fig. D.1.
In the text

Fig. D.3. Best-fit SEDs and P(z)s for U-53105 (top), A-26130 (middle) and A-76468 (bottom). The top section of each plot shows 1.5” × 1.5” NIRCam cutouts for the source. The left inset displays the observed photometry (black circles) with 1σ upper-limit uncertainties (black triangles). The best-fit SED derived with BAGPIPES is shown as a red curve with a 0.1 dex confidence interval (shaded red area). The right inset reports the normalized redshift probability distributions from different SED-fitting codes: the P(z) of Bagpipes in red; the CIGALE posterior in yellow; and the mean EAZY posterior averaged over all templates described in Sect. 5 (purple). The vertical dashed black line marks the z = 8 threshold.

In the text

	Fig. D.3. Continued.
In the text

Fig. E.1. Signal in the convolved long-wavelength NIRCam cutouts is shown for our stacked UHR candidate sample, with insets corresponding to F277W (top left), F356W (top right), F410M (bottom left), and F444W (bottom right). The blue line represents the signal, in μJy, for each pixel in the stacked image as a function of its distance from the center of the cutout (measured in pixels). A dashed green line depicts the Gaussian fit used to estimate the aperture photometry radius, while a dashed red line indicates the mean background level. The vertical dashed blue lines mark the inferred radius for fiducial aperture photometry.

In the text

Fig. G.1. Effective radii (in kiloparsecs) vs. redshift for our UHR candidate sample in the F277W, F356W, and F444W bands, respectively. We report their physical sizes estimated with galight vs. their z > 8 best-fit redshift solutions as colored stars. As a comparison, we report the empirical effective radius vs. redshift scaling retrieved in Westcott et al. (2025) as a continuous black line with 1σ (gray) and 2σ (light gray) confidence intervals.

In the text

	Fig. G.2. Same as Fig. G.1, but for our three HELM candidates and CURION.
In the text