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

The first year of JWST observations has revealed a wealth of surprises, including the remarkable overabundance of luminous galaxies in the epoch of reionization (EoR) relative to earlier expectations (Bouwens et al. 2015;Finkelstein 2016;Stark 2016;Finkelstein et al. 2022a;Robertson 2022). The Hubble Space Telescope (HST)ʼs discoveries at z > 8 told a story of a Universe rapidly growing at z ∼ 9 (Oesch et al. 2018), yet very few candidates were identified at z ∼ 10. This suggested that galaxies grow in lock step with their halos with roughly equal starforming efficiencies at all times, where "efficiency" is here defined as the effective stellar baryon fraction, ò å = M å / ( f b M halo ), where M å is the stellar mass, f b = 0.156 is the cosmic baryon fraction (Planck Collaboration et al. 2020), and M halo is the halo mass. Within the first few hundred million years (at z > 10) the dearth of galaxy candidates from the pre-JWST era was considered to be a natural consequence of the halo-growthlimited conversion of baryons into stars (Bagley et al. 2024;Harikane et al. 2023), a process that was thought to be independent of redshift at early times (e.g., Tacchella et al. 2013;Mashian et al. 2016;Stefanon et al. 2017;Oesch et al. 2018;Bouwens et al. 2023; though some work has suggested its evolution; e.g., Coe et al. 2013;McLeod et al. 2015;Finkelstein 2016;McLeod et al. 2016;Finkelstein et al. 2022b).

Nevertheless, HSTʼs deepest surveys led to the discovery of GN-z11 (Oesch et al. 2014;Skelton et al. 2014;Oesch et al. 2016), then a z ∼ 11 candidate that was not only the most distant galaxy candidate identified before JWST's launch, but also one of the most luminous, with an observed rest-frame UV magnitude of M UV ∼ -21.5. Selected from 0.2 deg 2 of aggregate deep HST imaging, its implied volume density was rather rare but difficult to constrain with a single source.

We now know GN-z11 to be at z = 10.60 thanks to JWST NIRSpec observations (Bunker et al. 2023) with a star formation rate (SFR) ∼ 20 M e yr -1 and stellar mass ∼ 10 9 M e (Tacchella et al. 2023)-all well in place within the first 400 Myr after the Big Bang. Beyond its extraordinary brightness, further JWST observations of GN-z11 reveal even more surprises: it exhibits Lyα in emission, has a possible damping wing of the intergalactic medium (IGM) observed as a Lorentzian absorption profile in Lyα (previously only seen as a signature of neutral IGM absorption in quasars; Miralda-Escudé 1998), as well as signatures of a candidate accreting supermassive black hole, Lyα halo, and possible nearby companions (Scholtz et al. 2023). These observations place new constraints on the early assembly of the highest density peaks in the cosmic web. Now that JWST is discovering new galaxies beyond z > 8 by the dozens (if not hundreds; Adams et al. 2023aAdams et al. , 2023b;;Finkelstein et al. 2023a;Franco et al. 2023;Harikane et al. 2023), we can directly assess whether or not GN-z11 is so unique (Mason et al. 2023). Most new discoveries have covered much fainter intrinsic luminosities with deep NIRCam data obtained over somewhat narrow fields of view (<100 arcmin 2 ), but the COSMOS-Web survey (GO#1727; Casey et al. 2023) is uniquely suited to the discovery of bright, rare sources in the EoR. In this paper we report the discovery of several extremely bright candidate galaxies beyond z  10 found in the first 0.28 deg 2 of NIRCam imaging data from COSMOS-Web. Though found in a similar survey area as the 0.2 deg 2 of HST imaging used to find GN-z11, HST could not have selected our candidates, as their detection relies on the extraordinary depth provided by JWSTʼs long wavelength (LW) imaging. Throughout, we use a Planck cosmology (Planck Collaboration et al. 2020), AB magnitudes (Oke & Gunn 1983), and a Chabrier initial mass function (IMF; Chabrier 2003).

We select this sample of z  10 candidates from the COSMOS-Web survey (GO #1727, PIs: J. Kartaltepe and C. Casey;Casey et al. 2023), a 255 hr imaging program covering a contiguous 0.54 deg 2 in four NIRCam filters (F115W, F150W, F277W, and F444W); in parallel, observations with MIRI in one filter (F770W) are obtained. We refer the reader to Casey et al. (2023) for a detailed description of the survey design. In this paper, we explore the first two epochs of COSMOS-Web data taken in January 2023 and April 2023, respectively. In total, the area surveyed in these two epochs is 0.28 deg 2 for NIRCam and 0.07 deg 2 for MIRI (covering 25% of the NIRCam mosaic).

COSMOS-Web NIRCam data reduction was performed using the JWST Calibration Pipeline (Bushouse et al. 2023) version 1.10.0, with the addition of several custom modifications also implemented in other works (e.g., Bagley et al. 2024). This includes subtraction of 1/f noise and the background. The Calibration Reference Data System pmap-1075 was used, corresponding to NIRCam instrument mapping imap-0252. Astrometry is anchored to Gaia EDR3 and bootstrapped from HST/Advanced Camera for Surveys F814W imaging (Koekemoer et al. 2007) and COSMOS2020 catalogs (Weaver et al. 2022). The normalized median absolute deviation of the astrometry is less than 12 mas for all filters. Mosaics with a 30 mas pixel scale are produced in stage 3 of the pipeline. A forthcoming paper (M. Franco et al. 2024, in preparation) will provide a more complete description of the COSMOS-Web NIRCam image processing. Similarly, we reduce MIRI data using the same pipeline with similar custom modifications; a more complete description of COSMOS-Web MIRI imaging (S. Harish et al. 2024, in preparation) will follow.

Beyond JWST, we use the wealth of multiwavelength data in the COSMOS field to vet detections presented in this paper, from the HST/F814W imaging from the original COSMOS Survey (Koekemoer et al. 2007;Scoville et al. 2007), the Spitzer COSMOS Survey (Sanders et al. 2007), Subaru Telescope Hyper Suprime-Cam (HSC) imaging (Aihara et al. 2022), and Ultra-VISTA imaging (McCracken et al. 2012), updated to the most recent release DR5. A full description of these data sets is provided in Weaver et al. (2022) and Casey et al. (2023). 30 We note that none of the targets included in our analysis are detected in the COSMOS2020 photometric catalog (Weaver et al. 2022), which used a deep CHI_MEAN detection (Drlica-Wagner et al. 2018) image constructed using UltraVISTA YJHKs plus HSC iz bands to extract source photometry. While some of the sources have a marginal signal in UltraVISTA imaging (detailed later in Table 1), they are not of sufficient signal-tonoise ratio (S/N) to have satisfied the COSMOS2020 analysis criteria. The z > 7.5 sources identified in COSMOS2020 were limited to those presented in Kauffmann et al. (2022), which have similar rest-frame UV luminosities (M UV < -21) to those presented here, but are generally at lower redshifts 7 < z < 10 and identified over a wider area. Similarly, none of our sources have significant emission in the millimeter (four of 12 sources have some Atacama Large Millimeter/submillimeter Array (ALMA) coverage, while all are covered by deep single-dish data from SCUBA-2; Simpson et al. 2019;J. McKinney et al. 2024, in preparation), X-ray (Civano et al. 2016), or radio (Jarvis et al. 2016;Smolčić et al. 2017).

We note that the MIRI data cover only 0.07 deg 2 of the 0.28 deg 2 covered to date in COSMOS-Web (25%); of the 15 candidates we discuss in this paper (whose selection is outlined in the next section), only three are covered by MIRI 7.7 μm pointings. Of those, only one is detected: COS-z12-2 at 7.8σ (the measured flux density is given later in Section 4). We note that detection at this threshold, as well as nondetection, is not particularly constraining to the spectral energy distribution (SED) fits, though the constraints are included in our analysis. Had any candidates been detected at 7.7 μm at higher significance (100σ), it would suggest they are more likely to be lower-redshift contaminants.

Photometric catalogs for COSMOS-Web imaging were constructed using the model-based photometric package Sour-ceXtractor++ (SE++ Bertin et al. 2020;Kümmel et al. 2020). A detection image is constructed using a CHI_MEAN combination (Szalay et al. 1999;Drlica-Wagner et al. 2018) of all four NIRCam filters (F115W, F150W, F277W, and F444W). After source detection, an optimal Sérsic model is fit using each individual NIRCam band. Then extraction is carried out using the best-fit NIRCam-derived Sérsic model on the JWST imaging, HST imaging, and all ground-based data sets described in Weaver et al. (2022). Model-based photometry allows simultaneous photometric measurements to be made on images with different PSFs without degradation. This allows us to fold in constraints from the deep ground-based data without loss of resolution (thus photometric precision) in our space-based data.

Alongside the primary SE++ photometry, we use SE "classic" (Bertin & Arnouts 1996) to perform aperture photometry on PSFhomogenized images from HST and JWST. PSF homogenization is performed using empirically measured PSFs built with PSFEx (Bertin 2011) and the pypher package (Boucaud et al. 2016), which computes a homogenization kernel between two different PSFs. All images are PSF homogenized to the F444W image. The same detection image is used for SE "classic" as was for SE++. A number of tests were performed to check the consistency of the model-based and aperture-based photometry, all of which will be described in a forthcoming paper (M. Shuntov et al. 2024, in preparation). Here we briefly describe key adjustments that were made to the uncertainties in the SE++ catalog for the purposes of this work. Our initial catalog construction with SE++ could only be achieved with weight maps, which accounts for read noise but fails to account for Poisson noise from the sky or sources, thus SE ++ underestimated errors for faint or undetected sources. We address this by independently deriving a measure of the Poisson sky noise using thousands of randomly placed circular apertures of varying diameter. We then compared the Poisson sky noise estimates to the standard deviation of flux density measurements for an idealized Sérsic profile (convolved with the PSF) placed randomly throughout the mosaic and found 0 3 diameter (for HST and JWST imaging) and 1″ diameter apertures (for groundbased imaging) captured the idealized model uncertainty well. This measure of the Poisson sky noise was then added in quadrature with the SE++-derived uncertainty to generate the estimates of the source photometry given in Table 1. For the purposes of this work, we adopt the model-based photometry but also provide aperture-based measurements for reference to the reader.

All sources in the SE++ and SE "classic" catalogs are then fit with photometric redshifts using the EAzY SED-fitting tool (Brammer et al. 2008) with a combination of the default Flexible Stellar Population Synthesis templates (Conroy & Gunn 2010;specifically QSF 12 v3) and bluer templates optimized for selecting less dusty galaxies in the EoR from Larson et al. (2023), which provide a number of models with variable Lyα escape fractions. For our initial runs we adopt the reduced Lyα template set; we run EAzY a second time on candidate z > 10 galaxies using the template set without Lyα, and perceive no change in the resulting redshift probability density distributions (PDFs). We do not use a magnitude prior in our EAzY runs, but instead adopt the default flat redshift prior.

The initial identification of candidate z > 10 galaxies is performed using the following criteria on the SE++ modelbased photometric catalog:

1. S/N F814W < 2, S/N F115W < 2, S/N F277W 5, and S/N F444W 5, 2. best-fit photometric redshift from EAzY of z a 8.5, 3. the integral of the redshift probability distribution above redshift 6 is p(z > 6) 0.95 and above redshift 8 is p (z > 8) 0.80, and 4. the magnitude in F277W is [F277W] 27.5.

Of the ∼340,000 sources identified in the 0.28 deg 2 , 620 sources fulfill these criteria, of which 340 (55%) are identified as hot pixels and are thrown out using size threshold criteria (i.e., sources with an SE++ area smaller than 100 modeled pixels and a radius smaller than 0 01). All remaining 280 candidates are visually vetted. Because we did not directly select sources using JWST color cuts, many are somewhat red ([F277W] -[F444W] > 1) as well as spatially extended. As a result they are likely to be at lower redshifts than z ∼ 10 and thus removed; 86 sources remain as plausible z  10 candidates. Their distribution in observed [F277W] against photometric redshift is shown in Figure 1 including some other well-known high-z sources in the literature including GN-z11 (Bunker et al. 2023;Tacchella et al. 2023). From these 86, we closely scrutinize sources that are estimated to have absolute UV magnitudes brighter than M UV ∼ -21 and z 10 as shown in Figure 1. There are 15 clearly separated from the population of fainter sources found at similar redshifts in COSMOS-Web. Note that all 15 sources also pass the S/N criteria of our selection using SE "classic" aperture photometry. We limit our analysis to this subsample only in this paper, focusing on galaxies at the most extreme margins of luminosity and redshift. Note. The photometry for all of our z  10 candidates extracted using two different methods: Source Extractor (SE) "classic", which is aperture-based photometry measured in 0 30 diameter apertures on pointspread function (PSF)-homogenized imaging for space-based data only, and from SourceXtractor++ (SE++) single-Sérsic-model-based photometry measured on images in their native resolution.

We then apply more stringent fits to the 15 z  10 candidates in our data set. First, we repeat the photometric redshift fitting of all candidates using LePhare (Arnouts et al. 2002;Ilbert et al. 2006), Bagpipes (Carnall et al. 2018), and EAzY. We generate optimal fits after running each tool twice: once with a flat redshift prior from 0 < z < 20 and a second time with a flat prior from 0 < z < 7 to generate best-fit low-redshift template fits. All codes use the full flux density constraints and their uncertainties in all broadband filters (rather than upper limits).

Our EAzY runs are slight modifications of the initial tests used to select the sample: we allow a wider range of possible redshift solutions with finer redshift sampling, and shift the adopted template set from Larson et al. (2023) to those with no Lyα emission to account for a presumably neutral IGM at z  10.

For LePhare optimal fits, we follow the methodology of Kauffmann et al. (2022) and briefly summarize it here. Bruzual & Charlot (2003) templates are used spanning a range of star formation histories (SFHs; exponentially declining and delayed-tau) as in Ilbert et al. (2015), with two different dust attenuation curves (Calzetti et al. 2000;Arnouts et al. 2013). Emission line fluxes are added following Saito et al. (2020), allowing variation in line strength by a factor of 0.3 dex from expectation (Schaerer & de Barros 2009). Note that Lyα emission is included in these fits, though not at significantly high equivalent widths to impact the fits to the photometry of these bright sources.

For Bagpipes optimal fits, we implement a delayed-τ SFH as a fraction of the age of the Universefoot_3 at redshift z plus a recent, instantaneous burst lasting between 1 and 100 Myr. We use a Calzetti dust attenuation law (Calzetti et al. 2000) with stellar population models from Bruzual & Charlot (2003). We allow the absolute magnitudes of attenuation A V to span 0-3 to capture a reasonable range of attenuation, ranging from unobscured to values more heavy attenuated than are seen in the integrated light of typical submillimeter galaxies (e.g., da Cunha et al. 2015). We note that restricting to A V < 3 is required to prevent fits using rather unphysical models, for example extremely dust-attenuated dwarf galaxies with A V ∼ 6 at low redshift (z < 1). We later discuss why such low-mass, extreme attenuation models are inconsistent with the sample.

We also fit model libraries of brown dwarf SEDs from Morley et al. (2012) and Morley et al. (2014) to the photometry, which span temperatures 200-2000 K and a range of surface gravities; none of our sample are well fit by brown dwarf templates. As discussed in Section 3.2, all sources in our sample are also spatially resolved, further reinforcing the idea that it is unlikely any are brown dwarfs.

We then repeat a number of these fits on the photometry extracted from the SE "classic" 0 3 diameter apertures on all constraining space-based imaging (i.e., HST/F814W and the JWST/NIRCam filters). We determine that differences in the measured photometry largely do not impact the results, though some nuances of the differences (and their effects on the redshift probability distributions) are discussed again later in Section 5.1. We adopt the model-based photometry from SE++ as our fiducial photometry and find that all fits (with EAzY, LePhare, and Bagpipes) estimate <3% of the redshift probability distributions is at z < 7.

In this work, we adopt the posterior distributions of physical parameters from our Bagpipes tests to describe the characteristics of the sample, such as stellar mass, 100 Myr averaged SFRs, M UV , the rest-frame UV slope β, and the absolute magnitude of attenuation A V . The motivation for such a choice is further discussed in Section 3.2. In addition, Bagpipes is the only code that directly provides posterior distributions for all physical parameters, giving valuable insight into the covariances and the nature of potential contaminants. For sources with a significant fraction of their redshift solutions beyond z > 15 (those presented in Section 4.3), we enforce an artificial cap at z ∼ 15 for their physical characteristics, realizing solutions between 15 < z < 20 are far less likely than 13 < z < 15. We show all of the SED fits for the sample of 15 candidates in Figure 2.

For each source, we quantify a normalized χfoot_4 metric, which we calculate as c c = N n 2 2 bands , where N bands is the number of effective bands available to us that are most directly constraining for 10 < z < 14 galaxies. Here we adopt N bands = 7, taking the five space-based bands (HST F814W, and JWST F115W, F150W, F277W, and F444W) and we count two ground-based filters whose depths are most useful for this work: UltraVISTA H and Ks. A future work on the UV luminosity function (UVLF) from COSMOS-Web will explore a more careful, quantitative definition of c n 2 .

Note that c n 2 is not a reduced χ 2 , which would also account for the degrees of freedom in model fits. 32 That is a complex problem in SED fitting, as often models (used to generate templates in SED fitting) have more tunable parameters than Figure 1. The distribution of bright candidate z  10 galaxies we identify in the first 0.28 deg 2 of COSMOS-Web (gray points). The galaxies described in this paper (gray stars) are drawn from this sample, focusing on the particularly luminous subset (galaxies with initial estimates of M UV < -21). Note that the conversion from [F277W] to M UV as plotted is only approximate, as exact conversions depend on the rest-frame UV slope. The gray region marks the parameter space not explored in this work ([F277W] > 27.5). Well-known literature sources are shown in purple at their measured spectroscopic redshifts. The redshifts of GL-z10 and GL-z12 are given in accordance with their marginal detections of [O III] from ALMA (Bakx et al. 2023;Yoon et al. 2023) and their photometric uncertainties (Castellano et al. 2022).

galaxies have photometric detections, and reducing the searchable parameter space by making well-motivated physical choices differs greatly between SED-fitting tools.

We motivate such a normalized χ 2 to have a quantity that allows it to be directly compared between surveys that have a variety of different filters used to select their samples. For Cutouts and SEDs of all sources in our sample. Cutouts are 3″ × 3″ and include, from left to right: a stack of HSC griz, HST F814W, a stack of UltraVISTA Y and HSC y, JWST F115W, UltraVISTA J, JWST F150W, UltraVISTA H and Ks, and JWST F277W and F444W. The plotted SED shows model-based photometry and 2σ upper limits for photometric points below 2σ significance (gray). We overplot the best-fit high-redshift solutions from EAzY (red), LePhare (blue), and Bagpipes (green), and their corresponding redshift PDFs between 8 < z < 20 in the inset plot. In gray we show the best-fit redshift solutions forced to z < 7 from EAzY (dashed), LePhare (solid), and Bagpipes (dotted). Synthesized photometry is shown in open boxes for two SEDs only for clarity: the EAzY-based high-z solution (red) and the LePhare low-z solution (gray). Normalized c n 2 values (see Section 3.1) are given in each panel for both the high-z and low-z solutions.

example, galaxies selected in the CEERS survey (Finkelstein et al. 2023a) given the seven bands that are constraining for COSMOS-Web z > 10 candidates. Given the unique nuances of each redshiftfitting tool, we only directly compare c n 2 for low-z and high-z fits using the same fitting code (e.g., EAzY low z to EAzY high z, LePhare low z to LePhare high z, and Bagpipes low z to Bagpipes high z). Thus the criterion for a robust candidate is

For transparency in our selection, Figure 2 shows all 15 z > 10 candidates with M UV  -21, though some of them have viable low-redshift solutions as measured via c D n 2 . This includes COS-z12-4, which we have multiple reasons to believe is at low z (detailed in Section 4), as well as COS-z13-3 and COS-z14-2, which are only detected in the F277W and F444W bands, and are discussed further in Section 4.3. Note that sources detected in only two bands are naturally more difficult to select cleanly using a χ 2 metric as most data constraints are upper limits. Nevertheless, these sources are discussed individually in Section 4.3 and again in Section 5.1 regarding the probability that they are low-redshift interlopers.

We conducted a simple test to assess the accuracy of the photometric redshift estimates made using the suite of best-fit SEDs for the whole sample from EAzY, LePhare, and Bagpipes. These sets of templates are assumed to sample the breadth of realistic SEDs for z > 10 galaxies. We shifted the SEDs forward and backward in redshift, modeled synthetic photometry with the characteristic noise of our data, and remeasured the photometric redshifts. We found a systematic offset toward higher redshifts using EAzY, while Bagpipes exhibited no systematic offset. This is similar to the Bagpipes and EAzY systemic offsets seen in the first epoch of COSMOS-Web z ∼ 9 -11 galaxies from Franco et al. (2023). We adopt the physical parameters and redshifts measured by Bagpipes for the rest of this work.

We measure the sizes of the sample using the F277W band using both GALFIT (Peng et al. 2002(Peng et al. , 2010) ) and GALIGHT (Ding et al. 2020). The F277W band is chosen as it provides the highest S/N measurements for all galaxies in the sample with the best spatial resolution. GALFIT uses a least-squares fitting algorithm while GALIGHT uses a forward modeling approach to find the optimum Sérsic fit to a galaxy's twodimensional light profile, after accounting for the PSF. We use average PSF images generated from our 2023 April epoch of data measured using PSFEx (Bertin 2011). The size measurements were broadly consistent between the two fitting techniques and in all cases spatially resolved; we provide the half-light radius measurements from GALFIT in Table 2.

Here we provide more detailed descriptions and relevant information about each of the 15 candidate galaxies at z > 10. We note that the sample may be delineated into roughly three subsets: exceptionally luminous 10 < z < 12 candidates (M UV < -21.5), bright 10 < z < 12 candidates (-20.5 > M UV > -21.5), and z > 13 candidates with M UV < -20.5. These three samples contain five sources each. The z > 13 candidates are only detected in two JWST bands, F277W and F444W, and thus are not as well constrained in terms of their physical properties. Photometry is provided for all sources in Table 1 and the positions and derived physical properties are given in Table 2.

In the descriptions that follow, we quote the relative probability that a source is a low-redshift interloper; we draw these probabilities from the EAzY photometric redshift fits using the SE++ photometry, whose distributions are presented in Figure 3 and later discussed as an ensemble in Section 5.1. We chose EAzY redshift PDFs due to their simplicity and straightforward adoption of a flat redshift (and magnitude) prior. However, a general caveat of interpreting redshift PDFs is that the exact amplitude of the redshift peak is quite sensitive to the adopted template set, range of physical parameters governing the input SEDs, and to small differences in the adopted photometry. So while generally multiple independent codes find consist peaks in the distribution (also shown in the inset panels in Figure 2), the integral under each peak is somewhat uncertain.

We also address the presence of on-sky neighbors as possible physical associations, checking for consistency with low-redshift solutions to our candidates (particularly because close neighbors on the sky are statistically more likely to be physically associated; Kartaltepe et al. 2007;Shah et al. 2020). We also assess the possible contribution of gravitational lensing from neighbors and find it to be insignificant in all cases (it is most significant for COS-z10-1, where we demonstrate a lensing magnification no greater than μ ≈ 1.01).

This source is detected in five bands and exhibits a 2.5 mag drop between F150W and the 2σ upper limit in F115W. COS-z10-1 is formally detected in UltraVISTA H-and Ks-band imaging in addition to the three reddest NIRCam filters. Its redshift is fit to z phot = 10.3-10.4 using EAzY and LePhare and slightly lower at z phot = 9.7 from Bagpipes. COS-z10-1 has two neighbors within 1 5: one 0 47 to the southwest and another 1 21 to the south-southwest. The former has a photometric redshift of z phot = 0.97, and the latter z phot = 2.04, both of which are inconsistent with the forced lowredshift solutions for COS-z10-1 at z low = 2.4, indicating that physical association with the neighbors (adopting the lowredshift solution) is unlikely. We estimate the maximum lensing that could occur from these foreground sources as μ ≈ 1.01 given the estimated masses of the foreground objects of 10 8 M e . There is no significant low-redshift peak in the redshift PDF at low redshift for COS-z10-1, whose integrated probability of being at z < 7 is <0.1%.

This galaxy is detected >3σ at high S/N in F150W, F277W, and F444W with marginal ∼2σ signals in UltraVISTA H and Ks. It has a substantial factor of 5.2 ± 0.2 flux density drop (1.7 mag) in the F150W filter, and another factor of 2.6 ± 0.2 drop to the 2σ upper limit in F115W. COS-z12-1 does not have an abrupt spectral break, though its photometry can be explained well with redshift z = 12.0-12.5. Forced low-redshift solutions produce photometric redshift estimates at z low < 1 for EAzY and LePhare and at z low = 3.2 for Bagpipes. Solutions at z low < 1 would require extreme emission line strengths for a relatively low-mass galaxy (with a specific star formation rate (sSFR) of 10 -7 yr -1 ), thus are less likely than the z low = 3.2 Bagpipes solution, though all are significantly less well fit to the data than the high-redshift solutions (p(z < 7) ∼ 0.9%). A foreground neighbor, located 1 2 to the southeast, is fit with a photometric redshift of z phot = 2.55 but sufficiently distant to not contaminate COS-z12-1's photometry; its photometry and z phot is also distinct from COS-z12-1's lowredshift solutions to not suspect physical association. COS-z12-1 is the most intrinsically luminous z > 12 galaxy found in our sample with M UV =- -+ 22.19 0.17 0.10 .

Detected in four bands (F150W, UltraVISTA H, F277W, and F444W, with a marginal ∼2.3σ signal in UltraVISTA Ks), COS-z12-2 has a factor of 5.4 ± 0.4 drop in flux density (1.5 mag) between UltraVISTA H and F150W, which demarcates the candidate's Lyman break. Its photometric redshift estimates span z phot = 11.9-12.1 with possible low-redshift solutions z low = 1-3. The low-redshift solution that is most plausible is likely the strong-line emitter at z = 3.03 found with LePhare, though not quite as extreme a drop around λ = 1.5 μm would be expected compared to observations. We also tested the reliance of the photometric redshift on the low-S/N UltraVISTA data; after removing the UltraVISTA constraints, a z ∼ 12 solution is still favored for COS-z12-2, though with increased probability of a low-redshift solution (31%). The low-redshift probability is decreased with the H-band detection in particular. There are no close neighbors within 1 5 of COS-z12-2, rendering its photometry clean from contamination. There is a source 1 9 distant to the northwest that has a photometric redshift from COSMOS2020 of z phot = 3.68, which is closer to but still statistically distinct from COS-z12-2's low-redshift solution. We note that, depending on the exact tuning of template sets for EAzY or the adopted range of physical parameters used in Bagpipes, as much as 2% of the redshift PDF sits at z ∼ 3 (and it goes as high as 7.4% using only aperture photometry from space-based constraints). We note that COS-z12-2 is also detected with MIRI at 7.7 μm, the only galaxy in our sample to have such a detection; its flux density is S 7.7μm = 288 ± 37 nJy. Though the MIRI constraint is not plotted in Figure 2, this detection is consistent with an approximately flat spectrum (in F ν ) from the near-infrared and is included in the SED fits as an additional constraint. While nominally one might think it has a significant impact on the stellar mass estimate (in the rest-frame optical), this particular measurement has no significant impact due to its low S/N; it does reduce the uncertainty on the stellar mass slightly though. COS-z12-2 is the second brightest galaxy in this sample.

This source is detected in four bands (F150W, UltraVISTA Ks, F277W, and F444W, with marginal ∼2σ emission in UltraVISTA H) with a factor of 3.1 ± 0.5 drop from UltraVISTA H to F150W. COS-z12-3 has a redder rest-frame UV slope than the other galaxies in this sample, which introduces more possible degeneracies with a low-redshift origin to its photometry. Nevertheless, the strong break at λ ∼ 1.5 μm-a factor of 6-10 in flux density (2-2.5 mag)argues for a high-redshift solution, and its photometric redshift is consistently fit to z phot = 11.5-12.2. Without the Ultra-VISTA photometry, COS-z12-3 would be difficult to identify. For example, by removing the UltraVISTA constraints and refitting its redshift, COS-z12-3 defaults to z ∼ 3 solutions; the most critical band that places the source at z ∼ 12 is Ks, detected at 3.5σ. Fitting the photometric redshift using the HSC, HST, UltraVISTA Ks, and NIRCam data results in the z ∼ 12 solution with a ∼30% chance of a low-z solution; this is reduced to 2.4% when including the limits in the other UltraVISTA bands, in particular H band, where the low-z solution would demand a brighter flux density from the contribution of [O III] 5007 Å. COS-z12-3 has no neighbors within 3 0.

The inferred absolute UV magnitude of COS-z12-3 is = - .22 , it is the reddest of the bright sources in this paper; this combined with the high SFR leads naturally to a hypothesis that it may be detected (or detectable) at millimeter wavelengths. COS-z12-3 is one of the few sources covered by existing archival ALMA data at 2 mm from program 2021.1.00705.S (PI: O. Cooper); it is undetected with measured rms ∼ 0.08 mJy, which sets an approximate 2σ upper limit to the dust mass at z = 11.5 of 4 × 10 8 M e and L IR  6 × 10 11 L e (SFR  90 M e yr -1 ); both limits are not sufficiently constraining to be in conflict with the measured A V from the Bagpipes fit. Detecting millimeter continuum in Note. The positions are measured from the detection image used for the SE++ and SE "classic" catalogs. We provide three photometric redshifts for each source from Bagpipes, EAzY, and LePhare, but note most of the derived properties -including M UV , β UV , SFR 100Myr , and M å -are measured from the best-fit Bagpipes posterior distributions. R eff is measured from F277W imaging using GALFIT. We catagorize subsets of our sample as described in the text in Section 4, and include three sources which were removed for further analysis on suspicion they are low-z contaminants; if confirmed as high z, their properties may reflect what is given in this table.

such sources may require 2 mm observations with a sensitivity ∼ 0.01 mJy, though we stress that due to the negative K-correction, dust continuum observations of such sources do not confirm or refute a low-or high-redshift solution; spectroscopy is necessary.

The fifth exceptionally bright candidate we identify is COS-z12-4. The primary limitation in characterizing COS-z12-4 is the proximity of two neighbor galaxies whose emission is spatially confused in ground-based data. While nominally fit to a photometric redshift of z phot = 12.3-12.7 and = M UV --+ 21.90 0.15 0.15 , the SE++-model-based measurements for COS-z12-4 claim a 3σ detection in the HSC z band, and imaging from the original COSMOS SuprimeCam i-band imaging (Capak et al. 2007) may show a ∼1.5sigma detection; however, on close inspection, both may be positive noise fluctuations from the neighbors' emission. The lack of deep, optical constraints with high-resolution imaging deems the source less secure. The neighbors are 0 47 away to the west and 1 0 away to the southwest and are fit with photometric redshifts of z phot = 4.65 and 4.50, respectively. The former is also detected in COSMOS2020 with a consistent photometric redshift of z phot = 4.8. Indeed, the LePhare low-redshift photometric fit to COS-z12-4 is consistent with the neighbors' redshifts, at z low = 4.65. Though the best high-z fit still has a formally lower c n 2 than this low-z solution, the consistent photometry with a low-z neighbor is sufficient to cast significant doubt on the high-z solution. The situation is similar to the environmental coincidence that argued CEERS-93316 was z ∼ 4.9 and not z ∼ 16 from Naidu et al. (2022a), later confirmed by Arrabal Haro et al. (2023). We therefore choose to remove COS-z12-4 from our analysis for the rest of this paper. Nevertheless, we provide its measured physical characteristics if it were at z > 10 in Table 2. We emphasize that obtaining a spectrum of COS-z12-4 is important in case it is indeed an ultrahigh-redshift source; in that case we will have learned a valuable lesson about the claimed completeness of our survey, the true occurrence rates of chance projections with low-z sources, and the abundance of ultrabright z ∼ 12 sources.

COS-z10-2 is one of the more intrinsically red galaxies in our sample and has a slightly bimodal morphology in the NIRCam LW bands. It is formally detected in three NIRCam bands. The photometric redshift estimates fall between z = 9.1 -10.1. The derived rest-frame UV slope for COS-z10-2 is β = -1.22 - + 0.16 0.22 , the second reddest of the sample; despite its red color, it passes all c D n 2 criteria thanks to the strength of its Lyman break, which is a factor of 6 in flux (2 mag). The source has no neighbors, but its double-component morphology may cast some doubt on the reliability of single-Sérsic-profile-based measured photometry. However, the aperture photometry for this source is fully consistent with the model-based results (albeit with more uncertainty in the redshift PDF, see Figure 3). COS-z10-2 has the highest derived stellar mass estimate of any Figure 3. A comparison of redshift PDFs fit using EAzY to three sets of photometry for each high-z galaxy candidate. The results for the SE++ model-based photometry ("M"), which includes both space-and ground-based photometry, are shown in dark red. Light orange shows the results for the SE "classic" photometry extracted in 0 3 diameter apertures from PSF-homogenized HST and JWST NIRCam images only ("AS" for aperture, space-based only). Dark orange shows the distributions after joining together the aperture-based photometry from SE "classic" with the model-based photometric measurements from the ground-based data ("ASG," aperture with space and ground). This last fit is meant to highlight the relative value and importance of the ground-based constraints, particularly those from deep Subaru/HSC imaging in the optical and UltraVISTA imaging in the near-infrared. Inset are the percentages of the redshift PDFs that lie at z < 7. In the majority of cases, the model-based photometry and aperture-based photometry produce similar redshift PDFs and the addition of the ground-based constraints to the aperturebased photometry dramatically reduces the integral of the redshift PDF below z < 7 for most sources.

in our sample with

10 M e (a consequence of its redder color in the rest-frame UV).

COS-z10-3 is detected in three NIRCam bands. COS-z10-3 has a close neighbor 0 78 to the southeast that has a photometric redshift of z phot = 2.32 that is also found in COSMOS2020 with a similar photometric redshift. We note that this is somewhat consistent with one of the three lowredshift solutions for COS-z10-3 found by EAzY of z = 2.37, which raises the probability this is a low-redshift interloper. However, we note that such a solution is a relatively poor fit compared to its corresponding high-redshift solution ( c D = 0.9 n 2 ), sufficiently high to remain in the sample.

COS-z10-4 is similarly detected in three NIRCam bands with a compact core and diffuse extended emission. There is no indication that it is detected in UltraVISTA imaging. The stack of HSC griz imaging displays a puzzling excess emission to the southwest; because this is not spatially coincident with the source's position within 1″, it is not of significant concern. COS-z10-4 has no close neighbors, and has redshift solutions spanning z phot = 10.2-11.2, with significantly better high-z fits than forced low-z solutions.

COS-z11-1 has the highest redshift of the "bright" subset and is detected in the three NIRCam bands only, with detections in F150W, F277W, and F444W. There is no evidence from the ground-based cutouts of significant emission above the background noise. It has no neighbors and is fit to redshifts z phot = 11.2-11.7. Using model-derived photometry, 0.4% of the redshift PDF is at z < 7, though a significantly higher z < 7 probability is found while using aperture-based photometry alone (62%). However, when adding the groundbased constraints to the aperture-based photometry, the z < 7 solutions occupy a lower percentage of the distribution at 26.3%. 1.48 0.39 0.27 . The low-redshift solutions occupy ∼0.4% of the redshift probability distribution, though we note that adopting the aperture-based photometry results in a larger fraction of the redshift PDF at z < 7: 7.9% when using aperture-based photometry alone, and 13.8% when using aperture-based photometry added with model-derived ground-based constraints (as seen in Figure 3). The addition of the ground-based constraints perhaps makes the distinction between low-redshift and high-redshift solutions difficult for COS-z11-2 because of the relative offset in flux between F150W and UltraVISTA H band. Despite the additional ambiguity surrounding COS-z11-2, we keep it in the high-z sample for further analysis.

In an effort to explore the most extreme subset of new discoveries within the EoR, we have also identified a sample of candidate galaxies at z phot > 13 from our existing imaging. As a natural consequence of the design of our survey, sources with z phot > 12.5 are ostensibly only detected in F277W and F444W. None are sufficiently bright to be detected with Spitzer at 3.6 μm (a wavelength not covered with JWST imaging in COSMOS-Web), UltraVISTA filters, or MIRI at 7.7 μm. Beyond the limited filter set, the wavelength gap between F277W and F150W presents another challenge, making it difficult to quote photometric redshifts more precise than Δz ≈ 1. Another consequence of the limited photometry is that the difference in c n 2 between low-and high-redshift fits is diminished: c n 2 is overall lower because only two bands have S/N > 3. Indeed, all sources in this category fail at least one c D n 2 cut (primarily EAzY and LePhare).

Our approach toward candidates in this regime is thus somewhat conservative. Of an initial set of 31 sources in our initial EAzY catalog fit to photometric redshifts z phot > 12.5 with [F277W] < 27.5, we down-select to five viable candidates in this paper. Sources are rejected from the sample because they fail all c D > 0.6 n 2 criteria (for EAzY, LePhare, and Bagpipes), have >5% probabilities for z < 7 solutions, they have either diffuse morphologies with radii R eff > 0 5, or they are spatially unresolved in F277W (R eff < 0 15), or they are sufficiently red ([F277W] -[F444W] < 0), all of which cast significant doubt on a z > 13 solution. As a natural consequence of this approach, the five remaining sources-COS-z14-2, COS-z13-1, COS-z13-3, COS-z13-2, and COS-z14-1span somewhat bluer colors than the parent population with -0.7 < [F277W] -[F444W] < -0.3. All have integrated probabilities of being at z < 7 less than 0.4% using the fiducial model-based photometry.

We note that COS-z13-2 and COS-z14-1 both have MIRI coverage at 7.7 μm though neither is detected. The measured flux densities we obtain for them using SE++-model-based photometry are S 7.7μm = 26 ± 31 nJy and S 7.7μm = 49 ± 31 nJy, respectively.

Of the five candidates, only two show significantly elevated low-redshift probabilities when using the SE "classic" aperture photometry: COS-z14-2 has 11.9% of its PDF at z < 7 and COS-z14-1 has 7.2% of its PDF at z < 7. When adding the aperture photometry with ground-based constraints the z < 7 probability is reduced in COS-z14-1 to 1.8%, but is much higher (98.4%) for COS-z14-2. Though this result is inconsistent with our model-based constraints on COS-z14-2, we remove it from the analysis in the discussion. We also remove source COS-z13-3 from further analysis because it fails all c D n 2 criteria, even though the integral of its z < 7 probability is <0.5%.

For the purposes of the discussion and ensuing physical characteristics, we only retain three sources in the z > 13 sample: COS-z13-1, COS-z13-2, and COS-z14-1, but we provide descriptions of all five. We continue to stress that this sample is overall less robust than the sources described in Sections 4.1 and 4.2 and all, including those we have thrown out on suspicion they are low z, require spectroscopic confirmation.

Here we present a more detailed discussion of the ramifications of these discoveries. First we present a discussion regarding low-redshift interlopers, a direct measurement of their volume density and contribution to the UVLF, and a summary of their measured physical characteristics. We then present a more detailed discussion of their stellar mass estimates and their implied star formation efficiencies within the Lambda cold dark matter (ΛCDM) paradigm. We then discuss the potential gas content of the systems and last raise the possibility of future observations constraining the host galaxies' halo masses, which could also inform constraints on alternate cosmological frameworks.

JWST observations of high-redshift galaxies so far have made clear that galaxies at z < 7 with strong emission line equivalent widths in the rest-frame optical can masquerade as ultrahigh-redshift (z > 10) galaxy candidates in JWST's broadband filters (Fujimoto et al. 2023;Naidu et al. 2022a;Zavala et al. 2023;McKinney et al. 2023a). While discussion of this phenomenon arose prominently via model fitting to CEERS-93316, an exceptionally exciting and bright z ∼ 16 candidate (Naidu et al. 2022a;Zavala et al. 2023;Donnan et al. 2023), recent NIRSpec follow up confirming a low-redshift solution (z ∼ 4.9) highlights the complexity and difficulty in selecting robust, ultrahigh-redshift candidates (Arrabal Haro et al. 2023). This type of strong emission line contaminant may only be a substantial concern for z ∼ 4.6-4.8 contaminants in fields where coverage exists in a majority of NIRCam broadband filters. At those redshifts, strong emission line sources may appear as Lyman-break candidates at z > 15. However, COSMOS-Web has fewer filters and so emission line contaminants are possible over a broader range of redshifts. Examples of best-fit low-redshift solutions (restricted to z < 7) are shown in gray in Figure 2.

Figure 3 shows three different derivations of the full redshift PDFs from 0 < z < 20 for each source measured using EAzY assuming flat redshift priors. We compare our fiducial modelbased photometric constraints to aperture photometry. In all cases, the probability of a z < 7 solution using the model-based photometry is <3% for each source. The source with the highest probability at z < 7 is COS-z12-3 with 2.4%, likely caused by its comparatively red color; the rest of the sample has P(z < 7) < 1%. We emphasize again that the SE++-modelbased and SE "classic" aperture-based photometry are independently measured: the first from images in their native resolution, and the second from PSF-homogenized images. Because the aperture photometry only includes measurements from five bands (HST F814W and the JWST NIRCam bands), its derived redshift PDFs are broader and generally show an increased probability for a low-redshift solution. This increased uncertainty in the redshift PDF can be attributed to the limited number of bands used to constrain the redshift. To demonstrate the importance of the ground-based photometry for the photometric redshifts, we construct a hybrid photometric catalog which marries the HST + JWST constraints of the aperture photometry to the model-derived ground-based observations. We note that for these z > 10 candidates, the ground-based constraints are almost all modestly constraining nondetections (or a handful of low-S/N detections in UltraVISTA). As a result, this hybrid catalog can be thought of as a sanity check on our model-based photometry, as we have replaced the most constraining, high-S/N bands (i.e., JWST NIRCam) with aperture-photometry measurements. In all but two cases, COS-z10-4 and COS-z14-2, those additional ground-based constraints significantly reduce the probability of a viable low-z solution over the space-based-only PDFs. COS-z10-4 is kept in the sample given that it formally passes our c n 2 criteria for inclusion and COS-z14-2 is removed as previously discussed in Section 4.3.

An important caveat of the previous paragraph is the adoption of flat redshift priors, which is a somewhat standard, though potentially flawed, literature convention. The on-sky surface density of galaxies as a function of redshift declines steeply with increasing redshift, regardless of brightness. For example, the on-sky surface density of a z ∼ 4 27th magnitude (AB) source is ∼300-500 times higher than similarly bright z ∼ 10 sources according to some of the latest compilations of the UVLF (e.g., Finkelstein 2016; Harikane et al. 2023). By this argument, any source with a low-redshift peak exceeding 0.3% may have as much of a 50% chance of being a true lowredshift source, and those with 3% may be 10 times more likely to be low z than high z. Nevertheless, such a thought experiment does not adequately account for variations in the intrinsic SEDs, the lack of detection in deep optical stacks, or the measured sizes of the sample. A dedicated study focused on best practices in photometric redshift fitting to such samples would be timely, though beyond the scope of this work.

Given the broader range of potential contaminants in COSMOS-Web, we explored if the ensuing derived parameters for such low-redshift solutions in M å -SFR-A V space were physical. Digging into the posterior distributions of the physical properties found from Bagpipes low-redshift fits, we find the redshift range of probable contaminants spans 1 < z < 4 with median redshift 〈z〉 = 3, stellar mass 〈M å 〉 = 10 9 M e , 〈SFR〉 = 10 M e yr -1 , and attenuation 〈A V 〉 = 1.7. The attenuation (thus reddening of the stellar continuum), combined with high SFRs, is what is needed to reproduce the photometry of a z > 10 Lyman break. Note that allowing A V to vary by up to 6 generates another cluster of possible solutions at z < 1 with A V ∼ 5; these are largely unphysical, as they would imply extreme attenuation and high SFRs in low-mass dwarf galaxies (Bisigello et al. 2023). If such a phenomenon were common, submillimeter number counts would likely be dominated by such sources, and they are not (Casey et al. 2014;Fujimoto et al. 2016).

Even with more reasonable limits set on A V < 3, we note that one may expect a nonnegligible fraction (20%-30%) of contaminants to be detectable in existing (sub)millimeter imaging in the field; we estimate this fraction by inferring A UV from A V (where A UV ≈ 2.6A V ), converting to IRX ≡ L IR /L UV , and thus scaling M UV to L IR . Sources above ∼10 12 L e would be detectable in existing SCUBA-2 and AzTEC maps in the field (Aretxaga et al. 2011;Casey et al. 2014;Simpson et al. 2019). None of our sample are detected above 3σ detection limits in those data sets.

The solid angle covered by COSMOS-Web to date is 0.28 deg 2 , implying a survey volume between 9.5 < z < 12.5 ∼4.5 × 10 6 Mpc 3 . This redshift interval brackets the more confident candidates that are detected in more than two bands and the average redshift for galaxies in the bin is near z ∼ 11. While our z > 13 candidates have redshift PDFs extending out to z ∼ 20, solutions at z > 15 are unlikely, so we cap the volume estimate relevant for those sources at ∼2.4 × 10 6 Mpc 3 , corresponding to 13 < z < 15 (the z ∼ 14 sample).

Figure 4 shows the direct contribution of detected sources (calculated using the 1/V max method) presented in this paper to the UVLF at z ∼ 11, with measurements also provided in Table 3. We calculate the contribution of these sources to the UVLF through 500 Monte Carlo draws from the posterior distributions of measured M UV and redshift values from Bagpipes, where sources may be counted in different magnitude bins in different draws or fall outside of the designated redshift range. These measurements have not been corrected for incompleteness, yet they do address potential contamination. Broadly, we find that the volume density of very bright galaxies found in COSMOS-Web is well aligned with expectation from other recent JWST-measured LFs at z > 10. Figure 4 shows measurements at z ∼ 10-12 from Donnan et al. (2023), Franco et al. (2023), Harikane et al. (2023), Leung et al. (2023), McLeod et al. (2024), and data from the CEERS collaboration (S. Finkelstein, private communication). The two functional double power-law fits shown are from Donnan et al. (2023) and Leung et al. (2023).

Most noticeable from our data is the relatively flat slope of the UVLF at the bright end; this could be caused by incompleteness in our lower-luminosity bin, which we do not correct in this work. Considering that the magnitude bin at M UV = -22 is relatively complete, our data disfavor a Schechter function fit to the UVLF similar to other work. A more thorough observational derivation of the UVLF from COSMOS-Web will follow in a forthcoming paper (M. Franco et al. 2024, in preparation) and it will include completeness simulations and an extended measurement down to the threshold detection limit of the survey.

We show in Figure 5 the distribution of sources presented in this paper in absolute UV magnitude, rest-frame UV slope (β), and stellar mass. We compare to other reported candidates in the literature summarized recently in Franco et al. (2023). The four brightest sources discussed in Section 4.1 are shown in green in each panel. Those four sources have luminosities well matched to and exceeding GN-z11 at similar redshifts. Given the wide area covered by COSMOS-Web to date, it is clear that we are sensitive to the discovery of more intrinsically luminous sources brighter than M UV ≈ -20.5 beyond z ∼ 10. At a fixed redshift, most JWST discoveries from deeper but narrower fields are about 1-3 mag (or 2.5-15 times) fainter than the sample presented herein.

The rest-frame UV slopes of this sample are a bit redder than most Lyman-break galaxies (LBGs) at this epoch (the median presented in the literature at z > 8 is 〈β〉 = -2.3 ± 0.5 and the median of this sample is 〈β〉 = -1.7 ± 0.5); one source, COS-z12-3, is a significant outlier with β = -0.6 while all others are bluer than β = -1.2. Our sample is likely more red than the other samples for a few reasons: as exceptionally bright sources, they tend to have higher estimated stellar masses (see the earlier measurements of this relation from Finkelstein et al. 2012;Tacchella et al. 2022). Those works measured a direct correlation between β and M å but not β and M UV (though the latter relationship is derived in Topping et al. 2023). With higher masses, it is more likely that the stellar population is generally older or the SFH more complex, such that there are either proportionally fewer O stars contributing to the restframe UV flux or the dust reservoirs in such early galaxies has built up enough to redden the UV a small amount (see Ferrara et al. 2023;Ziparo et al. 2023). While more enhanced metallicity may also account for relatively red β slopes compared to low-metallicity galaxies with similar SFHs, a β above -2.0 points to either dust attenuation or less recent star formation as the cause of the flatter UV slope. We do issue some caution in the interpretation of the M å -β relationship, as both quantities are derived from SED fitting and have nonnegligible covariance.

While bluer UV slopes are found in our highest-redshift (z > 13) sample, we caution that this is likely a selection effect: galaxies with redder values of β at z ∼ 13 would have significant degeneracies with low-redshift solutions and thus (Adams et al. 2023a;Bouwens et al. 2023;Donnan et al. 2023;Finkelstein et al. 2023b;Franco et al. 2023;Harikane et al. 2023;Leung et al. 2023;McLeod et al. 2024). The two functional fits shown are the Donnan et al. (2023) z ∼ 10 double power-law fit (lighter blue) and the Leung et al. (2023) z ∼ 11 double power-law fit (darker blue). The tabulated UVLF data from this paper are provided in Table 3 along with a measurement at z = 14 (not shown here). The z ∼ 10 measurement from the first epoch of COSMOS-Web data in Franco et al. (2023) is shown in light green. These measurements do not account for incompleteness, but do account for contamination by adjusting the contribution of each source by the probability that it is indeed z  10. Our measurements are in agreement with other literature estimates of the UVLF at similar redshifts. Note. The measured contribution of our candidates to the UVLF at z = 11 and z = 14 have not been corrected for incompleteness.

be removed from our sample for failing our selection criteria. This bias, which is prevalent in the selection of LBGs at all redshifts, cannot easily be addressed without significant investments in spectroscopy for large samples.

Figure 6 shows the sizes of galaxies in the sample against stellar mass surface density and SFR surface density. All sources are spatially resolved with average sizes 〈R eff 〉 ∼ 500 pc, which reduces concern that their emission may be dominated by an active galactic nucleus (AGN; though an unresolved morphology would not guarantee it). Stellar mass and SFR surface densities are calculated by dividing the total M å or SFR by 2 to account for the M å or SFR internal to R eff , then we divide by pR eff 2 . The stellar mass surface densities are very similar to some of the most compact local elliptical galaxies (Lauer et al. 2007;Hopkins et al. 2010), and ultracompact dwarfs and superstar clusters (Haşegan et al. 2005;Evstigneeva et al. 2007;Hilker et al. 2007;McCrady & Graham 2007) though our sample is about 10 times larger in R eff . This may suggest that, as observed, these galaxies could be viable progenitors of elliptical galaxy cores with similar densities. The star formation surface densities are very similar to those seen in some z > 7 bright LBGs (e.g., Bowler et al. 2017), some local starbursts (e.g., the Great Observatories All-sky Luminous Infrared Galaxy Survey sample; Armus et al. 2009;McKinney et al. 2023b) and highredshift submillimeter galaxies (Hodge et al. 2016;Burnham et al. 2021), though the latter systems tend to be much larger (R eff > 1 kpc) with higher SFRs.

Mass derivations from rest-frame UV data are naturally uncertain. Nevertheless, JWST provides a longer wavelength lever arm than HST, into the rest-frame optical, to constrain the SEDs of galaxies even beyond z > 10. The few bands and lack of coverage in the rest-frame near-infrared miss the gold standard of stellar mass derivation. However, the young age of the Universe at z > 10 (<500 Myr) significantly narrows the dynamic range of possible SFHs of luminous LBGs. This places reasonable limits on the mass-to-light ratio and thus the underlying stellar mass of individual sources, with derived uncertainties ∼ 0.2 dex despite the lack of constraints in the Figure 5. The derived physical characteristics of our luminous z > 10 candidates relative to other samples of early Universe galaxies in the literature (gray points, the majority of which are from the JADES sample in the left panel; Hainline et al. 2024), explicitly highlighting the COSMOS-Web sample found by Franco et al. (2023) in dark gray. The four galaxies that are exceptionally luminous are shown as green stars, the bright 10 z 12 sample is shown as dark blue stars, and the z > 13 sample is shown as light blue stars; this color scheme persists in Figures 6, 7, and9. Left: the rest-frame absolute UV magnitude against redshift in comparison to the literature samples. We highlight three other luminous 10 < z < 12 sources: GN-z11 (now spectroscopically confirmed at z = 10.60; Oesch et al. 2016;Bunker et al. 2023), GL-z10, and GL-z12 (Naidu et al. 2022b;Castellano et al. 2022) with tentative spectroscopic identifications from ALMA (Bakx et al. 2023;Yoon et al. 2023). Right: the rest-frame UV slope β with stellar mass. Our sample is the reddest subset of candidates reported in the literature (Finkelstein et al. 2012;Tacchella et al. 2022;Topping et al. 2023). The dotted-dashed line is the best-fit relation derived for z = 8 galaxies from Finkelstein et al. (2012). No trend is detected in M UV vs. β (see Topping et al. 2023).

Figure 6. The F277W-measured half-light radii of the sample plotted against the stellar mass surface density (top panel) and star formation surface density (bottom panel). At the top, we overplot the average stellar mass surface densities of local elliptical galaxies from Lauer et al. (2007) and ultracompact dwarfs and superstar clusters from the compilation in Hopkins et al. (2010). In addition, we overplot the measured sizes and surface densities of the compact z ∼ 8 sources first discovered by Labbe et al. (2023) and whose sizes were analyzed in Baggen et al. (2023). The star formation surface densities are on par with other z > 7 bright LBGs (Bowler et al. 2017), local starbursts (Armus et al. 2009;McKinney et al. 2023b), and submillimeter galaxies with measured sizes (with R eff > 1 kpc; Burnham et al. 2021). The sources from our sample are colored by subsample as in Figure 5.

rest-frame optical. Given the potential implications of stellar masses exceeding 10 9 M e at z > 10 (every source in our sample exceeds this limit), we address possible sources of uncertainty in our stellar mass derivation here.

We first check the sensitivity of the adopted SFH on the resultant stellar mass; generally, star formation that occurs further in the past will have a higher mass-to-light ratio implied for the rest-frame UV, and thus a higher stellar mass. We compare our fiducial Bagpipes model, which superimposes a delayed-τ SFH with a recent, constant starburst to a model with only a delayed-τ SFH. The delayed-τ-only SFH models effectively increase the stellar mass estimates for fixed photometry by factors of 0.4-3 (see also Michałowski et al. 2012;Mitchell et al. 2013;Michałowski et al. 2014). Conversely, we can ask what fraction of stellar mass in our sample assembles during the recent, constant starburst phase; in the majority of cases, 99% of the stellar mass is attributed to a recent burst (forming within the previous 50 Myr on average). If we allow for an even more extreme and recent burst without the contribution from the delayed-τ model, the stellar masses are only reduced by ∼10% below the fiducial estimates, well within the reported uncertainties. In that sense, this illustrates that the stellar masses we derive in this work are a conservative lower limit for a normally behaved stellar population, one comprised of Population II and I stars with low, but not extremely low, metallicity.

If Population III stars dominate the light in the rest-frame UV, their expected top-heavy IMF (Hirano et al. 2014(Hirano et al. , 2015) ) would result in a UV continuum dominated by nebular, rather than stellar, emission. This would lead to a higher light-to-mass ratio from the UV (e.g., Schaerer & de Barros 2009) by factors ∼0.5-0.6 dex, reducing the highest mass estimates in our sample at z ∼ 12 from ∼5 × 10 9 M e to ∼10 9 M e . This would, of course, imply that first generation star formation dominates the energy output of these systems, which may be a difficult and somewhat extreme boundary condition, even for such highredshift galaxies given that their estimated halo masses are quite high, M halo ∼ 10 11 M e . Another likely consequence of the rest-frame UV being dominated by Population III stars would be a much bluer slope, β, than measured here.

Lastly, we consider the impact of AGN on our stellar mass estimates. It has become clear that actively accreting supermassive black holes may be a fairly prominent source of energy output at early times, and JWST has enabled the identification of AGN with lower-mass black holes out to earlier times (Larson et al. 2023). For AGN to impact the masses in our sample significantly, they would need to dominate the restframe UV luminosity by a factor of several over the stellar contribution. This would effectively translate to a lower limit on the AGN luminosity of M UV < -20.5 or L UV  1 -2 × 10 10 L e . At an Eddington ratio of ò ∼ 0.1, this would imply a minimum black hole mass of M BH ∼ 5-6 × 10 6 M e . Such masses would be somewhat unexpectedly large for the downward-revised stellar mass estimates of their host galaxies, 10 9 M e , so we assess this outcome as less likely than our stellar mass estimates without a significant AGN contribution. While some high-z supermassive black holes do seem unusually massive for their host galaxies (Kocevski et al. 2023;Larson et al. 2023), it does not appear that they contribute significant emission to the continuum.

The stellar masses of our sample, as estimated using our fiducial Bagpipes burst + delayed-τ model, are shown against redshift in Figure 7. On average their uncertainties are ∼0.2 dex; these are smaller than the average stellar mass estimates for other similar HST-selected or JWST-selected galaxies in the EoR. Again this is due to the age of the Universe at z > 10 being <500 Myr and the luminosities of these sources being so bright; these two restrictions severely limit the range of allowable SFHs (with conventional IMF assumptions), demanding steep growth via recent bursts, and with limited SFHs come more refined limits on stellar masses.

In stellar mass space, a subset of our sample pushes the bounds of the most massive sources that could plausibly be found in deep JWST surveys at z > 10; in particular, COS-z12-1, COS-z12-2, and COS-z12-3, with masses M å ∼ 4-10 × 10 9 M e at z ≈ 12. We overplot curves of constant number density in Figure 7 implied from the halo mass function (Sheth & Tormen 1999) scaled by the cosmic baryon fraction and a reasonable "efficiency" of ò å ≈ 0.1, where ò å represents the fraction of baryons that have been converted into stars within a halo over its integrated lifetime, ò å ≡ M å /( f b M halo ). We note that the typical peak of the stellar mass-to-halo mass relation (SMHR) is M å /M halo ∼ 1-3 × 10 -2 across a range of redshifts (Mandelbaum et al. 2006;Shankar et al. 2006;Conroy & Wechsler 2009;Behroozi et al. 2010Behroozi et al. , 2019;;Shuntov et al. 2022), implying a maximum ò å ∼ 0.2 (and one might expect this fraction to be lower at much earlier cosmic times where constraints on the SMHR do not yet exist). We also overplot Boylan-Kolchin (2023). The volume in the current COSMOS-Web data set is ∼few × 10 6 Mpc 3 in a bin of width Δz = 2. This comparison highlights the particularly unusual nature of the three massive galaxies at z ∼ 12, COS-z12-1, COS-z12-2, and COS-z12-3, whose existence defies expectation. Their presence demands either a higher abundance of massive halos at z ∼ 12 or an enhanced stellar baryon fraction ò å ∼ 0.2-0.5, implying that the star-forming efficiency is elevated (ò SF ∼ 1) for a significant fraction of the galaxies' SFHs. The possible SFHs of these systems are shown in Figure 9 and more details on the implications of their masses are explored in the discussion. The sources from our sample are colored by subsample as in Figure 5. curves of constant peak height ν, a measure of the fraction of mass contained in halos above a given mass threshold; higher values of ν mark increasingly massive halos where ν = 4.5 at z = 0 corresponds to a halo of mass ≈ 5 × 10 15 M e . This makes clear that the implied evolution of our sample is extreme with ν > 6 given ò å = 0.1: they represent the highest possible mass overdensities that will grow to host massive galaxy clusters in the present-day Universe.

In Figure 8, we directly calculate the cumulative stellar mass density in two bins centered at z = 10 and z = 12 with width Δz = 2. We use the full posterior distributions in z and M å as derived using Bagpipes for these calculations and assume Poisson uncertainties. Direct comparison to the halo mass function implies reasonable to high stellar baryon fractions with ò å ≈ 0.1-0.3 at z = 10, not too distinct from expectations from the SMHR (Behroozi et al. 2019). However, at z = 12, our discoveries imply higher stellar mass fractions, with ò å ∼ 0.2-0.5 (see also the forthcoming publications by K. Chworowsky et al. 2024, in preparation;M. Shuntov et al. 2024, in preparation). Harikane et al. (2023) present a comprehensive overview of early Universe discoveries from JWST in its first year, and they find that particularly bright galaxies (M UV < -19.5) found in deeper, smaller volume surveys require high efficiencies, 33 where

3, to produce stellar masses ∼ 10 8-9 M e at z = 12-16. Such candidates were not expected before JWST. Are such high stellar baryon fractions of ò å ∼ 0.2-0.5 realistic? Some theoretical work suggests they are. In particular, galaxies at z > 10, largely embedded in the neutral Universe before reionization, would not be bombarded with a background of UV radiation and thus the rapid collapse of molecular clouds could see very high rates of star formation (Susa & Umemura 2004). The feedback-free starburst (FFB) scenario presented in Dekel et al. (2023) provides a detailed account from first principles of how such a starburst might be powered; in such systems, the freefall time is <1 Myr and rapid star formation occurs before massive stars develop winds and supernova feedback occurs, and the external UV background has not yet been established. This is very similar to prior simulation work which has demonstrated certain regimes where feedback fails to regulate star formation (Torrey et al. 2017;Grudić et al. 2018).

With a delayed onset of feedback, one might expect the instantaneous star-forming efficiency (

) for short periods (5 Myr), to lead to to appreciably larger ò å of order a few tenths. Galaxies residing in halos of mass ∼ 10 11 M e at z ∼ 10 may be expected to have FFB gas densities, thus they might have stellar masses as high as ∼10 10 M e , SFRs in the tens of solar masses per year, and blue, compact (subkiloparsec) morphologies. This describes the properties of the z ∼ 12 massive subsample well: 〈M å 〉 ≈ 5 × 10 9 M e , 〈SFR〉 ≈ 40 M e yr -1 , and R eff ≈ 500 pc. Future spectroscopy of such targets may further clarify the applicability of the FFB model to such systems, particularly in measurements of metallicity and rest-frame UV slope (and thus presence of dust).

Another theoretical interpretation of the very luminous, early systems is provided in Ferrara et al. (2023), who suggest that short-lived bursts of super-Eddington star formation may blow out the majority of dust in early galaxies (with sSFR > 20 Gyr -1 ), making it possible to detect very blue, luminous galaxies beyond z > 10. The median estimated specific SFR in our sample calculated on a 100 Myr (10 Myr) timescale is á ñ = - + sSFR 10 100 4

9 Gyr -1 (á ñ = - + sSFR 15 10 8

12 Gyr -1 ). This contrasts with, e.g., the very blue but similarly luminous z ∼ 7 -8 candidates identified by Topping et al. (2022) with sSFR  30 Gyr -1 . Such systems may have relatively high stellar masses M å ∼ 10 8-9 M e , metallicity Z ≈ 0.1 Z e , and blue rest-frame UV slopes (β ∼ -2.6). Ziparo et al. (2023) describe that either dust ejection by radiation pressure could result in such blue slopes, or alternatively, a patchy interstellar medium with spatially distinct regions of obscured and unobscured stellar light. Our sample is not quite as blue (with median 〈β〉 = -1.7 ± 0.5) as their fiducial model, which demonstrates some inconsistency with the super-Eddington dust ejection hypothesis, but the patchy attenuation model may indeed be applicable to these systems. Future ALMA observations will provide crucial information as to the dust content of such systems.

Figure 9 shows another rendering of the stellar masses in our sample against redshift. Here we have directly shown what the progenitor population of the three most massive systems at z ∼ 12 might look like, via posteriors (inner 68% confidence interval) on the SFHs, or cumulative stellar mass growth. The stellar mass growth in these systems is overwhelmingly Figure 8. The cumulative stellar mass volume density as a function of stellar mass calculated at z = 10 and z = 12 from our sample. The halo mass curves, with three different integrated stellar efficiencies or stellar baryon fractions of ò å = 0.1, 0.32, and 1 are derived using the methodology outlined in Boylan-Kolchin (2023). The points are derived using the full posterior distributions in z and M å for each source in bins of width Δz = 1 centered on either redshift. At z = 10 we infer stellar baryon fractions ò å ∼ 0.1-0.3, while at z = 12 we infer higher stellar baryon fractions, ò å ∼ 0.2-0.5. Open circles represent the extrapolated sum of stars and molecular gas; the gas masses are inferred from the Kennicutt-Schmidt (KS) relation. At the highest masses at z = 12, it is apparent that all baryons could be comprised of stars and molecular gas, leaving little room for substantial reservoirs of, e.g., atomic gas. This highlights the need to gather cold gas observations of the sample. dominated by a recent burst, such that their stellar masses have grown at rates that far outpace the growth of their parent dark matter halos (at a fixed volume density). As pointed out in Section 5.4, such rapid and recent stellar mass growth provides the most conservative stellar mass estimates for the sample as a whole. The suggested burst-driven nature of COS-z12-1, COS-z12-2, and COS-z12-3 in particular would imply that their host galaxy halo masses may intrinsically be lower than one might expect given their stellar masses, on par with some of the less luminous galaxies in our sample. Recent work from the FIRE simulation (Sun et al. 2023) has shown that, indeed, no special adjustments are needed to reproduce the observed characteristics of very bright early JWST discoveries, which they find are driven by stochastic and recent bursts of star formation.

We should note that in most, if not all, theoretical models built to explain very massive z > 10 galaxies, rapid bursts of star formation happen on short timescales, 5 Myr. Our SED fitting using Bagpipes allows bursts to have either a very short or relatively long duration, up to 100 Myr. Unfortunately limitations in the current data set do not enable meaningful direct constraints to be placed on the burst timescale (i.e., the posterior distribution of the ages of the burst component is flat). However, it may be possible and indeed necessary that such exceptionally luminous systems have experienced a series of short bursts that are fairly well modeled by one long-duration burst lasting up to 100 Myr.

The ascent of such massive z ∼ 12 systems is so rapid that the population we identify at z > 13-with stellar masses an order of magnitude lower-could plausibly serve as a progenitor population, despite the short timescale (<70 Myr) between the two epochs. At later times, it is possible that galaxies like COS-z12-1, COS-z12-2, and COS-z12-3 evolve to become some of the Universe's first massive galaxies (Carnall et al. 2023;Glazebrook et al. 2023).

Most baryons in galaxy halos beyond z > 3 should be contained in gas and not stars (Walter et al. 2020). What ramifications do such high stellar masses in our sample have on the potential to observe their reservoirs of molecular and atomic gas? Such observations may yet prove crucial to our interpretation of their masses, and thus efficiencies.

First, it is worth recognizing that typical observed efficiencies in the star formation process rarely exceed ò SF ∼ 0.1 (Evans et al. 2009;Bigiel et al. 2010;Kennicutt & Evans 2012). In the context of the galaxies' gas supply, the star-forming efficiency is as ò SF ≡ SFR/M gas (the inverse of the gas depletion timescale) and normalized to 10 8 yr, thus represents the fraction of gas consumed every 100 Myr. This is not the same as ò å , which we call the stellar baryon fraction in this work, but others refer to it as the star formation efficiency; ò å may be thought of the integral form of ò SF . Similarly

b halo is also not the same as ò SF , as the later captures only baryonic processes. If we approximate M gas ≈ f b M halo (which is a firm upper limit to M gas ) with SFR ≈ 50 M e yr -1 , we find the average star-forming efficiency would need to be ò SF 0.32, in excess of limits seen in local Figure Stellar mass against redshift for the candidates identified in this paper (stars). Here we highlight the integrated SFHs for the three most massive z ∼ 12 candidates, which incorporate a delayed-τ plus recent starbursts: COS-z12-1 (dark green), COS-z12-2 (teal), and COS-z12-3 (purple). The dashed lines show stellar mass growth histories if the SFH is instead assumed to be delayed-τ only without a burst; this results in higher stellar masses that would have built up more mass early (z > 16). With the delayed-τ plus burst model, the stellar masses may have increased by an order of magnitude in less than 100 Myr, from z ∼ 14 progenitors that look much like the candidates we identify in Section 4.3. Other high-z candidate galaxies from the literature are shown in gray points. Stellar masses that are formally disallowed in ΛCDM are noted in dark gray, corresponding to the stellar mass threshold for the most massive halo in the full sky assuming ò å = 1. We also show the same threshold corresponding to 0.28 deg 2 , the area covered by COSMOS-Web in this work. Similarly, the most massive galaxy anticipated in all of COSMOS-Web, calculated by Lovell et al. (2023), is shown in the thin gray line; it assumes the halo baryon-to-stellar conversion efficiency varies, with average 〈ò å 〉 ∼ 0.06, and the confidence intervals about that mass threshold encompass the most massive z ∼ 12 galaxy (COS-z12-3) within ∼2σ.

molecular clouds but necessary to build the observed stellar populations.

Though untethered to direct observations at z > 10, we can alternatively estimate gas masses in our sample using the star formation surface density to gas mass surface density conversion, or the KS relation (Schmidt 1959;Kennicutt 1998). Adopting a unimodal KS relation with a power-law index of roughly 2 (Ostriker & Shetty 2011;Narayanan et al. 2012) and with measured star formation surface densities ranging 6 < Σ SFR < 100 M e yr -1 kpc -2 , we extrapolate that the molecular (H 2 ) gas masses would range from ∼2 × 10 8 -4 × 10 9 M e . This presumes that the size of the gas reservoir is similar to the stellar reservoir. This would imply molecular gas fractions (H 2 to total baryonic mass) ∼40% on average. The effect of the added component of molecular gas on the total baryonic mass volume density is shown in Figure 8 in open circles. In the case of the z ∼ 12 candidates, this shows that summing the stellar mass and molecular gas mass fully compensates for the predicted baryonic content of these early halos, leaving little room for other baryonic contributions, for example, like atomic hydrogen, H I, which is an essential building block of molecular gas and transitional state for primordial gas to be transformed into stars.

Follow-up ALMA observations of these systems will prove invaluable to provide an independent estimate of the galaxies' total baryonic budgets. For example, spectral scans for [O III] (at rest-frame 88 μm), will not only provide much needed spectroscopic confirmation of these sources but also facilitate a direct dynamical mass estimate, accounting for both gas and stars. In the case of efficient star formation, we might expect the dynamical mass constraint from [O III] to be approximately equal to f b M halo ≈ 1-2 × 10 10 M e . In the case of a lower efficiency and higher halo mass (as one may expect from alternative cosmological models, see the next section), one would expect the dynamical mass constraint to be factors of several times higher due to the overall larger baryonic mass present in the halo.

An alternative interpretation to the very high stellar baryon fractions implied in our sample (with ò å ≈ 0.2-0.5) is that the six-parameter ΛCDM model underestimates the number density of massive halos at early times. This revision to the cosmological framework could be explained through the early dark energy (EDE) model (Karwal & Kamionkowski 2016;Poulin et al. 2018), which suggests that an early episode of dark energy injection, near the time of matter-radiation equality (followed by ΛCDM evolution), could both explain the higher perceived abundance of massive halos at early times (Klypin et al. 2021) and resolve recent measurements of the Hubble constant tension (e.g., Riess et al. 2022). As discussed in Boylan-Kolchin (2023), the enhanced matter density, power spectrum slope, and σ 8 in EDE could even be used to explain the high stellar masses measured for some of the most massive candidates found to date by JWST (Labbe et al. 2023) over much smaller regions of the sky (though we also note that more recent work has suggested downward revision of their stellar mass estimates due to the contribution of strong emission lines and/or AGN; Endsley et al. 2023a;Labbe et al. 2023). Indeed, EDE predicts the most profound differences for the most massive halos, which we are sensitive to detecting in COSMOS-Web; at z  10, EDE predicts ∼10 times the number of halos >10 11 M e than expected by ΛCDM. Though not directly measurable in the data we present in this work, future follow-up spectroscopy may be able to place more meaningful constraints on the dynamical masses of these bright z > 10 sources, thus giving more direct measurements of the abundance of massive halos at early times. 5.8. Stochastic Bursts Driving M UV < -21 Galaxies in the Epoch of Reionization

Provided ΛCDM still holds, the most straightforward explanation for the presence of such extraordinarily massive galaxies at z > 10 is their rapid growth through stochastic bursts of star formation where a significant fraction of available baryons is efficiently cooled, condensed, and transformed to stars on <100 Myr timescales. The similar volume densities measured for very massive galaxies (∼10 10 M e ) relative to those 10 time less massive suggests that very high stellar baryon fractions (ò å ∼ 0.2-0.5) are not typical of the broader population; the steepness of the halo mass function would otherwise demand that sources 10 times less massive are ∼100 times more common. With burst-driven star formation, galaxies can deviate to higher ò å , and Malmquist bias ensures they are the first to be characterized.

This hypothesis is consistent with burst-driven star formation dominating the bright end of the UVLF as suggested in cosmological simulations (Shen et al. 2023;Sun et al. 2023). Such efficient and quick growth would be facilitated by the lack of a UV background in the prereionization era. These bursts may be comprised of several brief (<5 Myr) episodes of super-Eddington star formation (e.g., Ferrara et al. 2023), consistent with the FFB model (Dekel et al. 2023), though future dust, gas, and spectroscopic observations may provide crucial tests for such dust-poor, low-metallicity models. Similar burst-driven rapid growth has been suggested for other, similar populations of UV-luminous galaxies recently discovered by JWST (Endsley et al. 2023b;Dressler et al. 2023;Looser et al. 2023). Ideally, direct mass constraints on every luminous z > 10 candidate identified may better inform their SFHs. Another key technique that might be used to constrain the stochasticity of bright galaxies in the EoR is a clustering analysis, as proposed in Muñoz et al. (2023) where the M UV < -21 population bias may provide insights into their host halo masses.

If the bright end of the UVLF (e.g., M UV < -21) is dominated by stochastic bursts with intrinsic timescales of <100 Myr, then a natural consequence is an evolution in the shape of the UVLF around z ∼ 8 from a double power law (at z > 8) to a Schechter function (at z < 8). 34 This would be the result if the UVLF is calculated in approximately fixed-width bins; for example, at z  8, a bin of width Δz ≈ 1 corresponds to a timescale less than 100 Myr, meaning bursts with timescales ∼50 Myr will be observed with a duty cycle of order unity (>0.5), whereas at z ∼ 6 the duty cycle would be substantially lower (0.25). Conversely, a careful analysis of the epoch marking the transition between a double power law 34 We note that several works have argued the UVLF is intrinsically a double power law down to z ∼ 4, though such works primarily relate to M UV ∼ -23 sources whose luminosities are likely attributable to AGN (Stevans et al. 2018;Adams et al. 2020;Finkelstein & Bagley 2022;Harikane et al. 2022), and here we discuss the -20  M UV  -22 regime not thought to be AGN dominated.

and Schechter function may help us directly constrain the characteristic burst timescale of very bright (M UV < -21) galaxies without the need to invoke dust or feedback to suppress the bright end of the UVLF.

We have presented 15 intrinsically luminous candidate galaxies at 10  z  14 with estimated UV absolute magnitudes spanning -20.5 > M UV > -22; three are identified as probable low-z contaminants and the remaining 12 are separated in three subsets: exceptionally bright 10 < z < 12 galaxy candidates with M UV < -21.5, bright 10 < z < 12 galaxy candidates spanning -20.5 > M UV > -21.5, and z > 13 candidates with M UV < -20.5, which are only detected in F277W and F444W with more uncertain physical characteristics. These sources were identified in the first 0.28 deg 2 area covered by the COSMOS-Web survey (Casey et al. 2023); their detection is only made possible by the exquisite sensitivity of the JWST NIRCam LW channels.

The rest-frame UV luminosities are among the brightest sources ever identified at these redshifts, on average 1-3 mag brighter than other newly identified JWST z > 10 galaxy candidates. Their rest-frame UV colors are slightly redder as well (with 〈β〉 = -1.7 ± 0.5), perhaps hinting at more complex underlying SFHs or the existence of early dust reservoirs that redden the stellar continua. All sources are spatially resolved with average R eff ∼ 500 pc. Their stellar mass surface densities are on par with local elliptical galaxies and ultracompact dwarf galaxies. Their star formation surface densities are similar to other exceptionally luminous z > 7 LBG candidates as well as local luminous infrared galaxies. Their stellar masses span 10 8.5 -10 10 M e with volume densities roughly of order 10 -6 Mpc -3 .

Four of the 12 robust z > 10 candidates have UV luminosities similar to GN-z11 at similar or higher redshifts. Three of these four sources, COS-z12-1, COS-z12-2, and COS-z12-3, test the limits of early stellar mass assembly with M å ∼ 5 × 10 9 M e at z ∼ 12. Given their implied stellar mass densities ∼ 10 4 M e Mpc -3 at z ∼ 12, we infer that ∼20%-50% of the baryons in their halos have been converted into stars (ò å ∼ 0.2-0.5, where ò å = M å /( f b M halo )). This requires either a very high star formation efficiency from very early times (ò SF ∼ 1 at z  16) with stellar mass growth that far outpaces dark matter growth of the underlying halos, or alternatively, a higher abundance of high mass halos that might be possible in alternatives to ΛCDM. We favor the first explanation, of rapid burst-driven growth in the stellar reservoirs, making it possible to build ∼10 10 M e of stars in less than 100 Myr. Such stochastic episodes of star formation may be responsible for the underlying shape of the bright end of the UVLF; a double power law could simply arise at z > 8 after accounting for the duty cycle of stochastic bursts at the highest redshifts compared to z ∼ 6-8.

While we have made a best effort to present secure z > 10 candidates in this paper, follow-up spectroscopy is crucial to confirm the extraordinary nature of these candidates. The facilities best equipped for that follow-up work are JWST itself -which could give a direct spectrum of the rest-frame UV and optical-and ALMA, which could be used to constrain the cold gas content in and around their halos. These sources could represent the brightest galaxies JWST will find in any field at z > 9 (unless another large field-of-view survey like COSMOS-Web is conducted in the future), and thus they serve as an important laboratory for the formation and evolution of the first bright galaxies, including the search for Population III stars, the onset of metal and dust production, as well as direct constraints on the neutral gas fraction at very early times.