COSMOS-Web: Intrinsically Luminous z ≳ 10 Galaxy Candidates Test Early Stellar Mass Assembly

We report the discovery of 15 exceptionally luminous 10 ≲ z ≲ 14 candidate galaxies discovered in the first 0.28 deg2 of JWST/NIRCam imaging from the COSMOS-Web survey. These sources span rest-frame UV magnitudes of −20.5 > M UV > −22, and thus constitute the most intrinsically luminous z ≳ 10 candidates identified by JWST to date. Selected via NIRCam imaging, deep ground-based observations corroborate their detection and help significantly constrain their photometric redshifts. We analyze their spectral energy distributions using multiple open-source codes and evaluate the probability of low-redshift solutions; we conclude that 12/15 (80%) are likely genuine z ≳ 10 sources and 3/15 (20%) likely low-redshift contaminants. Three of our z ∼ 12 candidates push the limits of early stellar mass assembly: they have estimated stellar masses ∼ 5 × 109 M ⊙, implying an effective stellar baryon fraction of ϵ ⋆ ∼ 0.2−0.5, where ϵ ⋆ ≡ M ⋆/(f b M halo). The assembly of such stellar reservoirs is made possible due to rapid, burst-driven star formation on timescales < 100 Myr where the star formation rate may far outpace the growth of the underlying dark matter halos. This is supported by the similar volume densities inferred for M ⋆ ∼ 1010 M ⊙ galaxies relative to M ⋆ ∼ 109 M ⊙—both about 10−6 Mpc−3—implying they live in halos of comparable mass. At such high redshifts, the duty cycle for starbursts would be of order unity, which could cause the observed change in the shape of the UV luminosity function from a double power law to a Schechter function at z ≈ 8. Spectroscopic redshift confirmation and ensuing constraints of their masses will be critical to understand how, and if, such early massive galaxies push the limits of galaxy formation in the Lambda cold dark matter paradigm.

ple open-source codes and evaluate the probability of low-redshift solutions; we conclude that 12/15 (80%) are likely genuine z > ∼ 10 sources and 3/15 (20%) likely low-redshift contaminants.Three of our z ∼ 12 candidates push the limits of early stellar mass assembly: they have estimated stellar masses ∼ 5 × 10 9 M ⊙ , implying an effective stellar baryon fraction of ϵ ⋆ ∼ 0.2 − 0.5, where ϵ ⋆ ≡ M ⋆ /(f b M halo ).The assembly of such stellar reservoirs is made possible due to rapid, burst-driven star formation on timescales <100 Myr where the star-formation rate may far outpace the growth of the underlying dark matter halos.This is supported by the similar volume densities inferred for M ⋆ ∼ 10 10 M ⊙ galaxies relative to M ⋆ ∼ 10 9 M ⊙ -both about 10 −6 Mpc −3 -implying they live in halos of comparable mass.At such high redshifts, the duty cycle for starbursts would be of order unity, which could cause the observed change in the shape of the UVLF from a double powerlaw to Schechter at z ≈ 8. Spectroscopic redshift confirmation and ensuing constraints of their masses will be critical to understanding how, and if, such early massive galaxies push the limits of galaxy formation in ΛCDM.

INTRODUCTION
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 expectation (Bouwens et al. 2015;Finkelstein 2016;Finkelstein et al. 2022a;Stark 2016;Robertson 2022).Hubble'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 star-forming efficiencies at all times (where 'efficiency' is here defined as an effective stellar baryon fraction, ϵ ⋆ = M ⋆ /(f b M halo ), where M ⋆ is the stellar mass, f b = 0.156 the cosmic baryon fraction, Planck Collaboration et al. 2020, and M halo 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-growth-limited conversion of baryons into stars (Bagley et al. 2022;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. 2015McLeod et al. , 2016;;Finkelstein 2016;Finkelstein et al. 2022b).
Nevertheless, Hubble's deepest surveys led to the discovery of GN-z11 (Skelton et al. 2014;Oesch et al. 2014Oesch 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 Hubble imaging, its implied volume density was rather rare but difficult to constrain with a single source.* NSF Graduate Research Fellow † NASA Hubble Fellow 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 ∼20 M ⊙ yr −1 and stellar mass ∼10 9 M ⊙ (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 observed as a Lorentzian absorption profile in Lyα (and 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; Finkelstein et al. 2023;Harikane et al. 2023;Franco 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. 2022) 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 Hubble imaging used to find GN-z11, Hubble could not have selected our candidates, as their detection relies on the extraordinary depth provided by JWST's long wavelength 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).∼ 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 MUV < −21).Note that the conversion from [F277W] to MUV as plotted is only approximate, as exact conversions depend on rest-frame UV slope.The gray region marks 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 [Oiii] from ALMA (Bakx et al. 2023;Yoon et al. 2023) and photometric uncertainties (Castellano et al. 2022) .
We select this sample of z > ∼ 10 candidates from the COSMOS-Web Survey (GO #1727, PIs Kartaltepe & Casey;Casey et al. 2022), a 255 hour 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. (2022) 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. 2022) version 1.10.0, with the addition of several custom modifications also implemented in other works (e.g.Bagley et al. 2022).This includes subtraction of 1/f noise and background.The Calibration Reference Data System (CRDS) pmap-1075 was used, corresponding to NIRCam instrument mapping imap-0252.Astrometry is anchored to Gaia-EDR3 and bootstrapped from the Hubble/F814W imaging (Koekemoer et al. 2007) and COSMOS2020 catalogs (Weaver et al. 2022).The normalized median absolute deviation on 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.) 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.) will follow.
Beyond JWST, we use the wealth of multiwavelength data in the COSMOS field to vet detections presented in this paper, from the Hubble/F814W imaging from the original COSMOS Survey (Scoville et al. 2007;Koekemoer et al. 2007), the Spitzer COSMOS Survey (Sanders et al. 2007), Subaru Telescope Hyper Suprime-Cam (HSC) imaging (Aihara et al. 2022), and UltraVISTA imaging (McCracken et al. 2012), updated to the most recent release DR5.A full description of these datasets is provided in Weaver et al. (2022) and Casey et al. (2022) 1 .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 χ 2 detection image constructed using UltraVISTA Y JHKs plus HSC iz bands to extract source photometry; 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 herein, 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 (4 of 12 sources have some ALMA coverage, while all are covered by deep single-dish data from SCUBA-2; Simpson et al. 2019, McKinney et al., in prep), X-ray (Civano et al. 2016), or radio (Jarvis et al. 2016;Smolčić et al. 2017).
We note that MIRI data covers 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 3 are covered by MIRI 7.7µm pointings.Of those, only one is detected: COS-z12-2 at 7.8σ (the measured flux density given later in § 4).We note that detection at this threshold, as well as non-detection, is not particularly constraining to the 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.

PHOTOMETRY, SOURCE SELECTION, & MEASUREMENTS
Photometric catalogs for COSMOS-Web imaging were constructed using the model-based photometric package SourceXtractor++ (SE++; Bertin et al. 2020;Kümmel et al. 2020).A detection image is constructed using a χ 2 combination (Szalay et al. 1999) of all four NIRCam filters (F115W, F150W, F277W, and F444W), which is then used for reference to extract photometry from each individual band.Extraction is carried out on JWST imaging, Hubble imaging, and all ground-based datasets described in Weaver et al. (2022).Model-based photometry allows simultaneous photometric measurements to be made on images with different point spread functions (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 Source Extractor (SE) 'classic' (Bertin & Arnouts 1996) to perform aperture photometry on PSFhomogenized images from Hubble and JWST.PSF homogenization is performed using empirically-measured PSFs built with PSFEx (Bertin 2011) and the pypher package (Boucaud et al. 2016) that 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 model-based and aperture-based photometry, all of which will be described in a forthcoming paper (M.Shuntov & L. Paquereau et al.).Allowing full flexibility in the SE++ catalog construction leads to some underestimated errors for faint or undetected sources in certain bands; we address this by setting a floor to the noise estimates in each filter equal to the shot noise measured in random circular apertures of diameter 0. ′′ 3 (for Hubble and JWST imaging) and 1 ′′ (for ground-based imaging).For the purposes of this work, we adopt model-based photometry but also provide aperture-based measurements for reference to the reader.
All sources in the SE++ and 'classic' catalogs are then fit to photometric redshifts using the EAzY SED-fitting tool (Brammer et al. 2008) with a combination of the default Flexible Stellar Population Synthesis (FSPS) templates (Conroy & Gunn 2010, specifically QSF 12 v3) and bluer templates optimized for selecting less dusty galaxies in the EoR from Larson et al. (2022), which makes a number of models available 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 with-out Lyα, and perceive no change in the resulting redshift probability density distributions (PDFs).We do not use a magnitude prior in our EAzY runs, thus adopt the default flat redshift prior.
The initial identification of candidate z > 10 galaxies is performed using the following criteria on the SE++ model-based photometric catalog: Of ∼340,000 sources identified in the 0.28 deg 2 , 620 sources fulfill this criteria, of which 340 (55%) are identified as hot pixels and are thrown out using size threshold criteria (i.e.sources with SE++ area less than 100 modeled pixels and radius smaller than 0.01 arcsec).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] magnitude 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 fifteen clearly separated from the population of fainter sources found at similar redshifts in COSMOS-Web.Note that all 15 sources also pass the SNR criteria of our selection using SE 'classic' aperture photometry.We limit analysis to this subsample only in this paper, focusing on galaxies at the most extreme margins of luminosity and redshift.
We then apply more stringent fits to the fifteen z > ∼ 10 candidates in our dataset.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 second with a flat prior from 0 < z < 7 to generate best-fit low redshift template fits.All codes use 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. (2022) 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 here.Bruzual & Charlot (2003) templates are used spanning a range of star formation histories (exponentially declining and delayed) 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τ star formation history as a fraction of the age of the universe2 at redshift z plus a recent, instantaneous burst lasting between 1-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 heavily attenuated than seen in typical submillimeter galaxies ' integrated light (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 photometry to model libraries of brown dwarf SEDs from Morley et al. (2012) and Morley et al. (2014), spanning temperatures 200 K-2000 K and a range of surface gravities; none of our sample are well fit by brown dwarf templates.As discussed in § 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.Hubble/F814W and the JWST/NIRCam filters).We determine that differences in 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 § 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 star formation rates, M UV , the restframe UV slope β, and the absolute magnitude of attenuation A V .The motivation for such a choice is further discussed in § 3.2.In addition, Bagpipes is the only code that directly provides posterior distributions in all physical parameters, giving valuable insight on covariances and the nature of potential contaminants.For sources with significant fractions of their redshift solutions fit beyond z > 15 (those presented in § 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.

Quantifying Goodness of Fit
For each source, we quantify a normalized χ 2 metric, which we calculate as χ 2 n = χ 2 /N 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 (Hubble F814W, and JWST F115W, F150W, F277W, and F444W) and we count two ground-based filters whose depths, ∼27.5-28.2AB at 5σ (Casey et al. 2022)   [nJy] [nJy] [nJy] [nJy] [nJy] [nJy] [nJy] [nJy] [nJy] [nJy] [nJy] [nJy] [nJy] COS-z10-1 0. .Note that χ 2 n is not a reduced χ 2 , which would also account for the degrees of freedom in model fits3 That is a complex problem in SED fitting, as often models (used to generate templates in SED fitting) have more tunable parameters than 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 could be directly comparable between surveys that have a variety of different filters used to select their samples.For example, galaxies selected in the CEERS survey (Finkelstein et al. 2023) are constrained with more deep space-based photometric bands, and would naturally have a higher value of χ 2 , thus larger differential values ∆χ 2 between low and high redshift solutions.Normalizing by the number of constraining bands would make selection of sources in these different surveys directly comparable.For example, the oft-used ∆χ 2 > 4 criterion (e.g.Finkelstein et al. 2023;Hainline et al. 2023) here would translate to ∆χ 2 n > ∼ 0.6 given the 7 bands that are constraining for COSMOS-Web z > 10 candidates.Given the unique nuances of each redshift fitting tool, we only directly compare χ 2 n 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 criteria for a robust candidate is if ∆χ 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 ∆χ 2 n .This includes COS-z12-4, which we have multiple reasons to believe is at low-z (detailed in § 4), as well as COS-z13-3 and COS-z14-2, which are only detected in the F277W and F444W bands, discussed further in § 4.3.Note that sources detected in only two bands are naturally more difficult to cleanly select using a χ 2 metric as most data constraints are upper limits.Nevertheless, these sources are discussed individually in § 4.3 and again in § 5.1 regarding the probability that they are low-redshift interlopers.

Other measurements
We conducted a simple test to assess the accuracy of the photometric redshift estimates using the suite of best-fit SEDs for the whole sample from EAzY, LePhare and Bagpipes.This set of templates is assumed to sample the breadth of realistic SEDs for z > 10 galaxies.We shifted these SEDs forward and backward in redshift, modeled synthetic photometry with characteristic noise of our data, and remeasured photometric redshifts.We found a systematic offset towards 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 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 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 two dimensional light profile, after accounting for the PSF.We use average PSF images generated from our April 2023 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.

SOURCE DETAILS
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 in 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 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 § 5.1.We chose EAzY redshift PDFs due to their simplicity and straightforward adoption of a flat redshift (and magnitude) prior.However, a gen- eral caveat of interpreting redshift PDFs is that the exact amplitude of redshift peaks are quite sensitive to the adopted template set, range of physical parameters governing the input SEDs, and to small differences in adopted photometry.So while generally multiple independent codes find consist peaks in the distribution (also shown for 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).

Exceptionally Luminous 10
This source is detected in five bands and exhibits a 2.5 magnitude 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 low-redshift solutions for COS-z10-1 at z low = 2.4, indicating that physical association with the neighbors (adopting the low-redshift 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 < ∼ 10 8 M ⊙ .There is no significant low-redshift peak in the redshift probability density distribution at low redshift for COS-z10-1, whose integrated probability of being at z < 7 is <0.1%.

COS-z12-1
This galaxy is detected at high signal-to-noise in three bands (F277W, F444W and UltraVISTA Ks) with a substantial 5× flux density drop (1.7 magnitudes) in the F150W filter, and another factor of 4 drop (1.5 magnitudes) 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 a 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 (sSFR > ∼ 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 highredshift solutions (p(z < 7) ∼ 0.9 %).A foreground neighbor, located 1. ′′ 2 to the southeast, is fit to 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.10 −0.17 .

COS-z12-2
Detected in five bands (F150W, UltraVISTA H and Ks, F277W and F444W), COS-z12-2 has a factor 4 drop in flux density (1.5 magnitudes) between UltraV-ISTA 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.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 COS-MOS2020 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 probability density distribution sits at z ∼ 3 (and it goes as high as 7.4% using only aperture photometry on space-based constraints).This makes COS-z12-2 slightly less robust as a z > 10 candidate than others in this exceptionally luminous category.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; this detection is consistent with an approximately flat spectrum (in F ν ) from the near-IR.COS-z12-2 is the second brightest galaxy in this sample.

COS-z12-3
This source is detected in five bands (F150W, Ultra-VISTA H and Ks, F277W and F444W) with a strong factor of 5 drop from UltraVISTA H to F150W.COS-z12-3 has a redder rest-frame UV slope than 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 magnitudes)-argues for a high-redshift solution, and its photometric redshift is consistently fit to z phot = 11.5 − 12.2.Without the Ul-traVISTA photometry, COS-z12-3 would be difficult to identify.COS-z12-3 has no neighbors within 3. ′′ 0.
The inferred absolute UV magnitude for COS-z12-3 is M UV = −21.58+0.07 −0.31 , and its dust-corrected SFR of 54 ± 16 M ⊙ yr −1 is the highest of the sample.With a restframe UV slope of β = −0.60 +0.22  −0.33 and A V = 0.58 +0.22 −0.33 , 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) in the millimeter.COS-z12-3 is one of the few sources covered by existing archival ALMA data at 2mm from program (2021.1.00705.S, PI: O. Cooper); it is undetected with measured RMS of ∼0.08 mJy, which sets an approximate 2σ upper limit to the dust mass at z = 11.5 of 4×10 8 M ⊙ and L IR < ∼ 6 × 10 11 L ⊙ (SFR < ∼ 90 M ⊙ 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 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 does not confirm or refute a low-or high-redshift solution; spectroscopy is necessary.

COS-z12-4
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 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 SuprimeCam i-band may show a ∼1.5σ 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 to 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 χ 2 n 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 ultra high-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 ultra-bright z ∼ 12 sources.4.2.Bright 9 < z < 13 Candidates 4.2.1.COS-z10-2 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 with marginal detections in Ul-traVISTA J, H, and Ks.The marginal J-band detection combined with the deep upper limit from F115W pegs the redshift right around z ∼ 9 at the low redshift end of our sample.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.22  −0.16 , the second reddest of the sample; despite its red color, it passes all ∆χ 2 n criteria thanks to the strength of its Lymanbreak, which is a factor of 6 in flux (2 magnitudes).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 in our sample with M ⋆ = (1.1 +0.7 −0.3 ) × 10 10 M ⊙ (a consequence of its redder color in the rest-frame UV).

COS-z10-3
COS-z10-3 is detected in three NIRCam bands; though formally undetected (<3σ) in UltraVISTA, cutouts clearly show excess emission consistent with the 2σ upper limit in H-band.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 low-redshift 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 highredshift solution (∆χ 2 n = 0.9), sufficiently high to remain in the sample.

COS-z10-4
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 UltraV-ISTA imaging.The stack of HSC griz imaging displays a puzzling excess emission to the south-west; 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
COS-z11-1 is 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, adding the ground-based constraints in with the aperture-based photometry, the z < 7 solutions occupy a lower percentage of the distribution at 26.3%.

COS-z11-2
COS-z11-2 is detected in three NIRCam bands plus UltraVISTA H-band, such that a clear break is detected between UltraVISTA H and F150W.It lacks neighbors and has redshift solutions spanning z = 10.6 − 11.9.The low-redshift solutions for this source cluster around z = 2.1 − 2.6 and are only marginally worse (∆χ 2 n = 0.3 − 0.5) than the high redshift solutions, likely because it is both relatively faint and redder than other sources in the sample, with M UV = −20.77+0.23 −0.34 and β = −1.48+0.27  −0.39 .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% using aperture-based photometry alone, and 13.8% using aperture-based photometry added with modelderived ground-based constraints (as seen in Figure 3).The addition of ground-based constraints perhaps makes the distinction between low-redshift and high-redshift solutions difficult in 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.

Candidates at z > 13
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 χ 2 n between low-and high-redshift fits is diminished: χ 2 n is overall lower because only two bands have SNR > 3. Indeed, all sources in this category fail at least one ∆χ 2 n cut (primarily EAzY and LePhare).
Our approach towards 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 [F 277W ] < 27.5, we downselect to five viable candidates in this paper.Sources are rejected from the sample because they fail all ∆χ 2 n > 0.6 criteria (for EAzY, LePhare and Bagpipes) and have >5% probabilities of 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 ([F 277W ] − [F 444W ] < 0) to 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-1-span somewhat bluer colors than the parent population with −0.7 < [F 277W ] − [F 444W ] < −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 using SE 'classic' aperture photometry: COS-z14-2 has 11.9% of the PDF at z < 7 and COS-z14-1 has 7.2% 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 future analysis in the discussion.We also remove source COS-z13-3 from further analysis because it fails all ∆χ 2 n 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 sources described in § 4.1 and § 4.2 and all, including those we have thrown out on suspicion they are low-z, require spectroscopic confirmation.

DISCUSSION
Here we present a more detailed discussion of the ramification 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 ΛCDM.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.

The Possibility of Low Redshift Interlopers
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 ultra high redshift (z > 10) galaxy candidates in JWST's broadband filters (Zavala et al. 2022;McKinney et al. 2023a;Fujimoto et al. 2022;Naidu et al. 2022a).While discussion of this phenomenon arose prominently via model fitting to CEERS-93316, an exceptionally exciting and bright z ∼ 16 candidate (Donnan et al. 2023;Zavala et al. 2022;Naidu et al. 2022a), recent NIRSpec follow-up confirming a low redshift solution (z ∼ 4.9) highlights the complexity and difficulty in selecting robust, ultra high-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 probability density distributions from 0 < z < 20 for each source measured using EAzY assuming flat redshift priors.We compare our fiducial model-based photometric constraints to aperture photometry.In all cases, the probability of a z < 7 solution using modelbased 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++ model-based 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 (Hubble F814W and the JWST NIR-Cam 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 Hubble+JWST constraints of the aperture photometry to the model-derived groundbased observations.We note that for these z > 10 candidates, the ground-based constraints are almost all modestly constraining non-detections (or a handful of low SNR 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-SNR bands (i.e.JWST NIR-Cam) with aperture photometric 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 spacebased only PDFs.COS-z10-4 is kept in the sample given that it formally passes our χ 2 n criteria for inclusion and COS-z14-2 removed as previously discussed in § 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 27 th 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 low redshift source, and those with 3% may be 10× more likely 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 distribution of physical properties 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 ⊙ , ⟨SFR⟩ = 10 M ⊙ yr −1 , and attenuation ⟨A V ⟩ = 1.7.The attenuation (thus reddening of stellar continuum), combined with high star formation rates, is what is needed to reproduce the photometry of a z > 10 Lyman break.Note that allowing A V to vary 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, Narayanan, & Cooray 2014;Fujimoto et al. 2016).
Even with more reasonable limits set on A V < 3, we note that one may expect a non-negligible fraction (20-30%) of contaminants to be detectable in existing (sub)millimeter imaging in the field; we estimate this fraction by inferring A U V from A V (where A U V ≈ 2.6A V ), converting to IRX ≡ L IR /L UV , and thus scaling M UV to L IR .Sources above ∼ 10 12 L ⊙ would be detectable in existing SCUBA-2 and AzTEC maps in the field (Aretxaga et al. 2011;Casey et al. 2014; (green stars) against recent JWST literature measurements (Harikane et al. 2023;Donnan et al. 2023;McLeod et al. 2023;Leung et al. 2023;Franco et al. 2023, Finkelstein et al., in preparation).The two functional fits shown are the Donnan et al. ( 2023) z ∼ 10 double powerlaw fit (lighter blue) and the Leung et al. (2023) z ∼ 11 double powerlaw 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.Simpson et al. 2019).None of our sample are detected above 3σ detection limits in those datasets.

Volume Density and UVLF Contribution
The solid angle covered by COSMOS-Web to-date is 0.28 deg 2 , implying the survey volume between 9.5 < z < 12.5 of ∼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 UV luminosity function at z ∼ 11, with measurements also provided in Table 3.We calculate the contribution of these sources Note-The measured contribution of our candidates to the UV luminosity function at z = 11 and z = 14 have not been corrected for incompleteness.
to the UVLF through 500 Monte Carlo draws from the posterior distributions of measured M UV and redshift 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 luminosity functions at z > 10.Most noticeable from our data is the relative 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 disfavors 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 (Franco et al., in prep) and it will include completeness simulations and an extended measurement down to the threshold detection limit of the survey.

Physical Properties of Bright z > 10 Candidates
We show in Figure 5 the distribution of sources presented in this paper in absolute UV magnitude, restframe 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 § 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 todate, 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 fixed redshift, most JWST discoveries from deeper but narrower fields are about 1-3 magnitudes (or 2.5 − 15×) 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 other samples for a few reasons: as exceptionally bright sources, they tend to have higher estimated stellar masses (see earlier measurements of this relation from Finkelstein et al. 2012;Tacchella et al. 2022).Those works measured direct correlation between β and M ⋆ but not β and M UV .With higher masses, it is more likely that the stellar population is generally older or the star-formation history more complex, such that there are either proportionally fewer O stars contributing to the rest-frame UV flux or the dust reservoirs in such early galaxies has built up enough to redden the UV a small amount (cf.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 star-formation histories, β 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 non-negligible 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 be removed from our sample for failing our selection criteria.This bias, which is prevalent in 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 star formation rate surface density.All sources are spatially resolved with average sizes ⟨R eff ⟩ ∼ 500 pc, which reduces concern that their emission may be dominated by AGN (though an unresolved morphology would not guarantee it).Stellar mass and star-formation rate surface densities are calculated by dividing the total M ⋆ or SFR   2023)), 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 shown as dark blue stars, and the z > 13 sample as light blue stars; this color scheme persists in Figures 6, 7, and 9. Left: The rest-frame absolute UV magnitude against redshift in comparison to 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 (Castellano et al. 2022;Naidu et al. 2022b) with tentative spectroscopic identifications from ALMA (Bakx et al. 2023;Yoon et al. 2023).Right: The rest-frame UV slope β with redshift.Our sample is the reddest subset of candidates reported in the literature (Finkelstein et al. 2012;Tacchella et al. 2022).The dot-dashed line is the best-fit relation derived for z = 8 galaxies from Finkelstein et al. (2012).No trend is detected in MUV versus β.
by two to account for the M ⋆ or SFR internal to R eff , then we divide by πR 2 eff .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 Ultra Compact Dwarfs and Super Star 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 galaxies' cores with similar densities.The star formation surface densities are very similar to those seen in some z > 7 bright Lyman-break Galaxies (e.g.Bowler et al. 2017), some local starbursts (e.g. the Great Observatories All sky LIRG Survey, GOALS, sample, Armus et al. 2009;McKinney et al. 2023b) and high redshift submillimeter galaxies (Burnham et al. 2021;Hodge et al. 2016), though the latter systems tend to be much larger (R eff > 1 kpc) with higher star formation rates.

Stellar Mass Uncertainties
Mass derivations from rest-frame UV data are naturally uncertain.Nevertheless, JWST provides a longer wavelength lever arm than Hubble, into the rest-frame optical to constrain the SEDs of galaxies even beyond z > 10.While the few bands and lack of coverage in the rest-frame NIR misses the gold standard of stellar mass derivation, the young age of the Universe at z > 10 (<500 Myr) significantly narrows the dynamic range of possible star-formation histories of luminous LBGs, and thus places reasonable limits on the mass-to-light ratio, thus the underlying stellar mass of individual sources.Given the potential implications of stellar masses exceeding 10 9 M ⊙ at z > 10 (all of our sample exceeds these limits), we address possible sources of uncertainty in our stellar mass derivation here.
We first check the sensitivity of the adopted star formation history (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 restframe UV, thus 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 lowski et al. 2012Micha lowski et al. , 2014;;Mitchell et al. 2013).Conversely, we can ask what fraction of stellar mass in our sample assembles during the recent, constant starburst phase; in a 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 ex- treme and recent burst without the contribution from the delayed-τ model, the stellar masses are only reduced by ∼10% below the fiducial estimates, well within 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 restframe 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 of ∼0.5-0.6 dex, reducing the highest mass estimates in our sample at z ∼ 12 from ∼ 5 × 10 9 M ⊙ to ∼ 10 9 M ⊙ .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 high redshift galaxies given that their esti-mated halo masses are quite high, M halo ∼ 10 11 M ⊙ .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 significantly impact the mass in our sample, they would need to dominate the rest-frame UV luminosity by a factor of several over the stellar contribution.This would effectively translate to a lower limit on AGN luminosity of M UV < −20.5 or L UV > ∼ 1 − 2 × 10 10 L ⊙ .At an Eddington ratio of ϵ ∼ 0.1, this would imply a minimum black hole mass of M BH ∼ 5 − 6 × 10 6 M ⊙ .Such masses would be somewhat unexpectedly large for the downward-revised stellar mass estimates of their host galaxies, < ∼ 10 9 M ⊙ , so we assess this outcome as less likely than our stellar mass estimates without 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.

Massive Beacons: the assembly of the first megalithic halos?
In stellar mass space, a subset of our sample push 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 ⊙ at z ≈ 12.We overplot curves of constant number density on 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 (Shankar et al. 2006;Mandelbaum et al. 2006;Conroy & Wechsler 2009;Behroozi et al. 2010Behroozi et al. , 2019;;Shuntov et al. 2022), implying 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 curves of constant peak height ν, a measure of the fraction of mass contained in  Kolchin (2023).The volume in the current COSMOS-Web dataset is a ∼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 higher abundance of massive halos at z ∼ 12 or an enhanced stellar baryon fraction ϵ⋆ ∼ 0.2 − 0.5, implying that the star forming efficiency be elevated (ϵSF ∼ 1) for a significant fraction of the galaxies' star formation histories.The possible star formation histories of these systems is shown in Figure 9 and more details on the implications of their masses explored in the discussion.Sources from our sample are colored by subsample as in Figure 5.
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 ⊙ .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 expectation from the stellar-to-halo-mass relation (Behroozi .The cumulative stellar mass volume density as a function of stellar mass calculated at z = 10 and z = 12 from our sample; halo mass curves, with three different integrated stellar efficiencies or stellar baryon fractions: ϵ⋆ = 0.1, 0.32 and 1, are derived using the methodology outlined in Boylan-Kolchin (2023).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; gas masses are inferred from the Kennicutt-Schmidt 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.et al. 2019).However, at z = 12, our discoveries imply higher stellar mass fractions, with ϵ ⋆ ∼ 0.2 − 0.5 (see also forthcoming publications by K. Chworowsky et al. and M. Shuntov et al.).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 4 , where ϵ = Ṁ⋆ /(f b Ṁhalo ) ∼ 0.3, to produce stellar masses ∼ 10 8−9 M ⊙ 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) 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 free-fall time is <1 Myr and rapid star formation occurs before massive stars develop winds and supernovae feedback 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 (ϵ = Ṁ⋆ /(f b Ṁhalo ) ≈ 1) for short periods of time ( < ∼ 5 Myr), leading to appreciably larger ϵ ⋆ of order a few tenths.Galaxies residing in halos of mass ∼ 10 11 M ⊙ at z ∼ 10 may be expected to have FFB gas densities, thus they might have stellar masses as high as ∼ 10 10 M ⊙ , SFRs in the 10's of M ⊙ yr −1 , and blue, compact (sub-kpc) morphologies.This describes the properties of the z ∼ 12 massive subsample well: ⟨M ⋆ ⟩ ≈ 5 × 10 9 M ⊙ , ⟨SFR⟩ ≈ 40 M ⊙ yr −1 , and R eff ≈ 500 pc.Future spectroscopy of such targets may further clarify applicability of the FFB model to such systems, particularly in measurements of metallicity and rest-frame UV slope (thus presence of dust).
Another theoretical interpretation of the very luminous, early systems is provided in Ferrara et al. (2023), which suggests 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.Such systems may have relatively high stellar masses M ⋆ ∼ 10 8−9 M ⊙ , metallicity Z ≈ 0.1 Z ⊙ , and blue rest-frame UV slopes (β ∼ −2.6).Ziparo et al. (2023) describes that either dust ejection by radiation pressure could result in such blue slopes, or alternatively, a patchy ISM 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.The median estimated specific SFR in our sample calculated on a 100 Myr (10 Myr) timescale is ⟨sSFR 100 ⟩ = 10 +9 −4 Gyr −1 (⟨sSFR 10 ⟩ = 15 +12 −8 Gyr −1 ).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 star formation histories, or cumulative stellar mass growth.The stellar mass growth in these systems are overwhelmingly 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 § 5.4, such rapid and recent stellar mass growth provides the most conservative stellar mass estimates for the sample as a whole.The suggested burstdriven 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 dataset do not enable meaningful direct constraints on burst timescale (i.e. the posterior distribution in ages of the burst component is flat).However, it may be possible that such exceptionally luminous systems have experienced a series of short bursts that are fairly well modeled by one long-duration burst 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 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).

Not all baryons are stars
Most baryons in galaxy halos beyond z > 3 should be contained in gas and not stars (Walter et al. 2020).Figure 9. Stellar mass against redshift for the candidates identified in this paper (stars).Here we highlight the integrated star formation histories (SFH) 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 § 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 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σ.

Disallowed in Λ CDM
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, 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 galaxies' gas supply, star-forming efficiency is defined ϵ 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 a stellar baryon fraction in this work, but others refer to as star-formation efficiency; ϵ ⋆ may be thought of the integral form of ϵ SF .Similarly ϵ ≡ Ṁ⋆ /(f 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 ⊙ yr −1 , we find the average star-forming efficiency would need to be ϵ SF ≥ 0.32, in excess of limits seen in local 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 Kennicutt-Schmidt (KS) relation (Schmidt 1959;Kennicutt 1998).Adopting a unimodal KS relation with powerlaw index roughly = 2 (Ostriker & Shetty 2011;Narayanan et al. 2012) and with measured star-formation surface densities ranging 6 < Σ SFR < 100 M ⊙ yr −1 kpc −2 , we extrapolate that the molecular (H 2 ) gas masses would range from ∼ 2 × 10 8 − 4 × 10 9 M ⊙ .This presumes that the size of gas reservoir is similar to the stellar reservoir.This would imply molecular gas fractions (H 2 to total baryonic mass) of ∼ 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, HI, 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 [Oiii] (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 [Oiii] to be approximately equal to f b M halo ≈ 1−2×10 10 M ⊙ .In the case of 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.

Alternative Cosmologies predict more massive halos early
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 (hereafter ΛCDM) 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 matterradiation equality (followed by ΛCDM evolution), could both explain a higher perceived abundance of massive halos at early times (Klypin et al. 2021) and resolve recent measurements of the Hubble 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 in some of the most massive candidates found to-date by JWST (Labbe et al. 2022) 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 contri-bution of strong emission lines and/or AGN; Endsley et al. 2023;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× the number of halos > 10 11 M ⊙ 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 dynamical masses of these bright z > 10 sources, thus give more direct measurements of the abundance of massive halos at early times.

Stochastic Bursts Driving M UV < −21 Galaxies in the EoR
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 measured for very massive galaxies (∼10 10 M ⊙ ) relative to those 10× 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× less massive are ∼100× 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 starformation dominating the bright end of the UVLF as suggested in cosmological zoom-in simulations (Sun et al. 2023).Such efficient and quick growth would be facilitated by the lack of UV background in the pre-reionization 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 Feedback-Free Starbursts model (Dekel et al. 2023), though future dust, gas and spectroscopic observations may provide crucial tests for such dust-poor, lowmetallicity models.Ideally, direct mass constraints on every luminous z > 10 candidate identified may better inform their star-formation histories.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 for their host halo masses.
If the bright end of the UVLF (e.g.M U V < −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) 5 .This would be the result if a 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 powerlaw 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.

CONCLUSIONS
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. 2022); 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 magnitudes brighter than other newly identified JWST z > 10 galaxy candidates.Their restframe UV colors are slightly redder as well (with ⟨β⟩ = −1.7 ± 0.5), perhaps hinting at more complex underlying star-formation histories 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 ultra compact 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 5 We note that several works have argued the UVLF is intrinsically a double power law down to z ∼ 4, though such work primarily relates to M UV ∼ −23 sources whose luminosity is likely attributable to AGN (Stevans et al. 2018;Adams et al. 2020;Harikane et al. 2022;Finkelstein & Bagley 2022), and here we discuss the −20 > ∼ M UV > ∼ − 22 regime not thought to be AGN dominated.
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 ⊙ at z ∼ 12.Given their implied stellar mass densities ∼10 4 M ⊙ 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 ⊙ 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 powerlaw 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 followup 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.Support for this work was provided by NASA through grant JWST-GO-01727 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.CMC, HA, MF, JM, ASL and ORC acknowledges support from the National Science Foundation through grants AST-2307006, AST-2009577, and

Figure 1 .
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 MUV < −21).Note that the conversion from [F277W] to MUV as plotted is only approximate, as exact conversions depend on rest-frame UV slope.The gray region marks 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 [Oiii] from ALMA(Bakx et al. 2023;Yoon et al. 2023) and photometric uncertainties(Castellano et al. 2022) Figure 2 -continued.

Figure 3 .
Figure 3.A comparison of redshift probability density distributions fit using EAzY to three sets of photometry for each highz 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 Hubble 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 of 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 aperture-based photometry dramatically reduce the integral of the redshift PDF below z < 7 for most sources.

Figure 4 .
Figure4.The direct contribution of sources presented in this paper to the UV luminosity function at z ∼ 11 (green stars) against recent JWST literature measurements(Harikane et al. 2023;Donnan et al. 2023;McLeod et al. 2023;Leung et al. 2023; Franco et al. 2023, Finkelstein et al., in  preparation).The two functional fits shown are theDonnan et al. (2023) z ∼ 10 double powerlaw fit (lighter blue) and theLeung et al. (2023) z ∼ 11 double powerlaw fit (darker blue).The tabulated UVLF data from this paper are provided in Table3along with a measurement at z = 14 (not shown here).The z ∼ 10 measurement from the first epoch of COSMOS-Web data inFranco 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.

Figure 5 .
Figure5.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. (2023)), explicitly highlighting the COSMOS-Web sample found byFranco et al. (2023) in dark gray.The four galaxies that are exceptionally luminous are shown as green stars, the bright 10 ≤ z ≤ 12 sample shown as dark blue stars, and the z > 13 sample as light blue stars; this color scheme persists in Figures6, 7, and 9. Left: The rest-frame absolute UV magnitude against redshift in comparison to 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(Castellano et al. 2022;Naidu et al. 2022b) with tentative spectroscopic identifications from ALMA(Bakx et al. 2023;Yoon et al. 2023).Right: The rest-frame UV slope β with redshift.Our sample is the reddest subset of candidates reported in the literature(Finkelstein et al. 2012;Tacchella et al. 2022).The dot-dashed line is the best-fit relation derived for z = 8 galaxies fromFinkelstein et al. (2012).No trend is detected in MUV versus β.

Figure 6 .
Figure6.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 top, we overplot the average stellar mass surface densities of local elliptical galaxies fromLauer et al. (2007) and Ultra Compact Dwarfs and Super Star Clusters from the compilation inHopkins et al. (2010).In addition, we overplot the measured sizes and surface densities of the compact z ∼ 8 sources first discovered byLabbe et al. (2022) and whose sizes were analyzed inBaggen 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).Sources from our sample are colored by subsample as in

Figure 7 .
Figure7.Stellar mass versus redshift for our sample relative to curves of fixed volume density, peak height, and other candidate galaxies reported in the literature (light gray points).Volume density curves (purple) and peak height (orange) are generated from expectation of the halo mass function scaled by the cosmic baryon fraction f b and an efficiency to convert baryons to stars of ϵ⋆ = 0.1 as inBoylan- Kolchin (2023).The volume in the current COSMOS-Web dataset is a ∼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 higher abundance of massive halos at z ∼ 12 or an enhanced stellar baryon fraction ϵ⋆ ∼ 0.2 − 0.5, implying that the star forming efficiency be elevated (ϵSF ∼ 1) for a significant fraction of the galaxies' star formation histories.The possible star formation histories of these systems is shown in Figure9and more details on the implications of their masses explored in the discussion.Sources from our sample are colored by subsample as in Figure5.
Figure8.The cumulative stellar mass volume density as a function of stellar mass calculated at z = 10 and z = 12 from our sample; halo mass curves, with three different integrated stellar efficiencies or stellar baryon fractions: ϵ⋆ = 0.1, 0.32 and 1, are derived using the methodology outlined inBoylan-Kolchin (2023).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; gas masses are inferred from the Kennicutt-Schmidt 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.
the UT Austin College of Natural Sciences for support.CMC also acknowledges support from the Research Corporation for Science Advancement from a 2019 Cottrell Scholar Award sponsored by IF/THEN, an initiative of the Lyda Hill Philanthropies.
Center (DAWN) is funded by the Danish National Research Foundation (DNRF) under grant No. 140.This work was made possible thanks to the CANDIDE cluster at the Institut d'Astrophysique de Paris, which was funded through grants from the PNCG, CNES, DIM-ACAV, and the Cosmic Dawn Center; CANDIDE is maintained by S. Rouberol.The French contingent of the COSMOS team is partly supported by the Centre National d'Etudes Spatiales (CNES).OI acknowledges the funding of the French Agence Nationale de la Recherche for the project iM-AGE (grant ANR-22-CE31-0007).MBK acknowledges support from NSF CAREER award AST-1752913, NSF grants AST-1910346 and AST-2108962, NASA grant 80NSSC22K0827, and HST-AR-15809, HST-GO-15658, HST-GO-15901, HST-GO-15902, HST-AR-16159, HST-GO-16226, HST-GO-16686, HST-AR-17028, and HST-AR-17043 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555.SG acknowledges financial support from the Villum Young Investigator grant 37440 and 13160.
on observations collected at the European Southern Observatory under ESO programmes 179.A-2005 and 198.A-2003 and on data obtained from the ESO Science Archive Facility with DOI https://doi.org/10.18727/archive/52,and on data products produced by CALET and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium.

Table 2 .
Sample Characteristics

Table 3 .
UV Luminosity Function Constraints