An Atlas of Color-selected Quiescent Galaxies at z > 3 in Public JWST Fields

We present the results of a systematic search for candidate quiescent galaxies in the distant universe in 11 JWST fields with publicly available observations collected during the first 3 months of operations and covering an effective sky area of ∼145 arcmin2. We homogeneously reduce the new JWST data and combine them with existing observations from the Hubble Space Telescope. We select a robust sample of ∼80 candidate quiescent and quenching galaxies at 3 < z < 5 using two methods: (1) based on their rest-frame UVJ colors, and (2) a novel quantitative approach based on Gaussian mixture modeling of the near-UV − U, U − V, and V − J rest-frame color space, which is more sensitive to recently quenched objects. We measure comoving number densities of massive (M ⋆ ≥ 1010.6 M ⊙) quiescent galaxies consistent with previous estimates relying on ground-based observations, after homogenizing the results in the literature with our mass and redshift intervals. However, we find significant field-to-field variations of the number densities up to a factor of 2–3, highlighting the effect of cosmic variance and suggesting the presence of overdensities of red quiescent galaxies at z > 3, as could be expected for highly clustered massive systems. Importantly, JWST enables the robust identification of quenching/quiescent galaxy candidates at lower masses and higher redshifts than before, challenging standard formation scenarios. All data products, including the literature compilation, are made publicly available.


INTRODUCTION
Over the last few years, the existence of a population of quenched and quiescent galaxies (QGs) at redshifts z ∼ 3 − 4 (e.g., Fontana et al. 2009;Straatman et al. 2014;Spitler et al. 2014) has been finally corroborated by the long sought after spectroscopic confirmations (Glazebrook et al. 2017;Schreiber et al. 2018a,b;Tanaka et al. 2019;Valentino et al. 2020;Forrest et al. 2020a,b;D'Eugenio et al. 2020D'Eugenio et al. , 2021;;Kubo et al. 2021;Nanayakkara et al. 2022).The combination of spectra and deep photometry have allowed for a first assessment of the physical properties of the newly-found early QGs.These properties include suppressed and minimal residual star formation rates (SFR), also supported with long-wavelength observations (Santini et al. 2019(Santini et al. , 2021;;Suzuki et al. 2022); emission from active galactic nuclei potentially pointing at a co-evolution with or feedback from their central supermassive black holes (Marsan et al. 2015(Marsan et al. , 2017;;Ito et al. 2022;Kubo et al. 2022); stellar velocity dispersions (σ ) and dynamical masses (Tanaka et al. 2019;Saracco et al. 2020) with possible implications on their initial mass function (Esdaile et al. 2021;Forrest et al. 2022); very compact physical sizes and approximately spheroidal shapes (Kubo et al. 2018;Lustig et al. 2021); and evidence that their large-scale environment may perhaps be overdense (Kalita et al. 2021;Kubo et al. 2021;McConachie et al. 2022;Ito et al. 2023).
Particular attention has been given to the reconstruction of the history (formation, quenching, and subsequent passive evolution) of distant QGs.A rapid and intense burst of star formation -compatible with that of bright sub-millimeter galaxies with depletion timescales of τ 100 Myr -is thought to drive the early mass assembly of the most massive and rarest systems (Forrest et al. 2020b) as established for z ∼ 2 QGs (Cimatti et al. 2008;Toft et al. 2014;Akhshik et al. 2022).However, a more steady stellar mass assembly at paces typical of galaxies on the main sequence at z > 4 might explain the existence of at least a fraction of the first QGs, likely less massive (Valentino et al. 2020).In this case, the population of dust-obscured "Hubble-dark" or "optically-faint" sub-millimeter detected sources could represent a good pool of candidate progenitors (Wang et al. 2019;Williams et al. 2019;Barrufet et al. 2022;Nelson et al. 2022;Pérez-González et al. 2022).These results stem from various approaches and their inherent uncertainties, such as the modeling of star formation histories (SFHs) with different recipes -parametric or not (Schreiber et al. 2018b;Ciesla et al. 2016;Carnall et al. 2018Carnall et al. , 2019;;Leja et al. 2019a;Iyer et al. 2019, K. Gould et al. in preparation), matching comoving number densities of descendants and progenitors, also including "duty cycles" (i.e., there have to be at least as many star-forming predecessors as quiescent remnants accounting for the time window in which such progenitors are detectable, Toft et al. 2014;Valentino et al. 2020;Manning et al. 2022;Long et al. 2022) or clustering analyses (Wang et al. 2019).
Debate continues on the exact mechanisms causing the cessation of the star formation at z 3 − 4, as well as at other redshifts (Man & Belli 2018 for a compendium).However, at high redshift there is the significant advantage of observing such a young Universe that classical "slow" quenching processes operating on ≥ 1 − 2 Gyr timescales at low redshifts are disfavored (e.g., strangulation or gas exhaustion Schawinski et al. 2014;Peng et al. 2015).Moreover, aided by sample selections favoring high detection rates over completeness, the best characterized spectroscopically confirmed QGs tend to show signatures of recent quenching (∼ a few hundred Myr) as in "post-starburst" galaxies rather than being prototypical "red and dead" objects (Schreiber et al. 2018b;D'Eugenio et al. 2020;Forrest et al. 2020b;Lustig et al. 2021;Marsan et al. 2022;Gould et al. 2023), even if examples of older populations are available (Glazebrook et al. 2017, M. Tanaka et al. in preparation).The analysis of larger samples of galaxies during or right after quenching could eventually help us understand the physics behind this phenomenon in the early Universe.
The exploration of post-quenching evolution is also in its infancy.There are indications of a simultane- Note-NIRCam depths: expressed as 5σ within the 0. 5 apertures used for the photometric extraction in the area covered by F150W / F200W / F277W / F356W / F444W (Appendix A). a The area covered by the group members has been masked (Appendix A).
b Effective area accounting for the gravitational lensing effect at z ∼ 3 − 5 (Section 4).
c No F150W coverage.
ous passive aging of the stellar populations and a rapid size evolution, but only modest stellar mass increase via dry minor mergers (Tanaka et al. 2019), resembling the second act of the popular "two-phase" evolutionary scenario that explains how z ∼ 2 QGs change over time (e.g., De Lucia et al. 2006;Cimatti et al. 2008;Naab et al. 2009;Oser et al. 2010).From the point of view of stellar dynamics, the small sample of QGs with available velocity dispersions does not allow for drawing any strong conclusions about possible evolutionary paths at constant or time varying σ yet (Tanaka et al. 2019;Saracco et al. 2020;Esdaile et al. 2021;Forrest et al. 2022).
These first results already paint a rich picture of how the earliest QGs formed and quenched and indicate several promising research avenues to explore.However, they relied on the availability of deep near-infrared photometry and ground-based spectroscopy, which come with obvious limitations on the spatial resolution and wavelength coverage.So far, these prevented us from unambiguously confirming if QGs exist at z > 4 (see Merlin et al. 2019;Carnall et al. 2020;Mawatari et al. 2020 for possible candidates), clearly defining the first epoch of sustained galaxy quenching, and ascertaining the existence of low-mass systems potentially quenched by different processes.
JWST enables us to break this ceiling, looking farther and deeper to catch the earliest QGs spanning a vast range of stellar masses.The first months of observations kept this promise and already offer a spectacular novel view on early galaxy evolution in general.In this work, we aim to capitalize on publicly available JWST multiwavelength imaging in 11 fields to find and quantify the population of early QGs, pushing the limits in redshift and mass affecting ground-based surveys.This paper is the first of a series addressing several of the contentious scientific points mentioned above.Here we will focus on the JWST -based selection of a robust sample of photometric QG candidates and on the bare-bones comoving number density calculations, taking advantage of the coverage of a relatively large combined area of ∼ 145 arcmin 2 at z ∼ 3 − 5 and the scattered distribution on the sky to reduce the impact of cosmic variance.Counting galaxies is a basic test for models and simulations and, in the case of distant QGs, it has generated quite some discussion on the robustness of current theoretical recipes (e.g., Schreiber et al. 2018b;Merlin et al. 2019).Also, accurate number densities are key ingredients to try to establish an evolutionary connection among populations across redshifts, thus affecting our view of the history of assembly of the first QGs.
The data collection, homogeneous reduction, and modeling are presented in Section 2. Our JWST -based color selection is described in Section 3, followed by the results on number densities contextualized within the current research landscape in Section 4. Throughout the paper, we assume a ΛCDM cosmology with Ω m = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 .All magnitudes are expressed in the AB system.All the reduced data, selected samples, and physical properties discussed in this work are publicly available online 1 .

DATA
In the following sections, we present our reduction and analysis of the photometric data.A dedicated paper will describe all the details of this process (G.Brammer et al., in preparation).The approach is similar to that in Labbe et al. (2022) andBradley et al. (2022), here including the recently updated zeropoints.

Reduction
We homogeneously process the publicly available JWST imaging obtained with the NIRCam, NIRISS, and MIRI instruments in 11 fields targeted during the first three months of operations (Table 1).We retrieved the level-2 products and processed them with the Grizli pipeline (Brammer & Matharu 2021;Brammer et al. 2022).Particular care is given to the correction of NIRCam photometric zeropoints relative to jwst 0942.pmap,including detector variations 2 .The results are consistent with similar efforts by other groups (Boyer et al. 2022;Nardiello et al. 2022) and with the more recent jwst 0989.pmapcalibration data.Corrections and masking to reduce the effect of cosmic rays and stray light are also implemented (see Bradley et al. 2022).For the PRIMER data, we introduce an additional procedure that alleviates the detrimental effects of the diagonal striping seen in some exposures.Finally, our mosaics include the updated sky flats for all NIR-Cam filters.We further incorporate the available optical and near-infrared data available in the Complete Hubble Archive for Galaxy Evolution (CHArGE, Kokorev et al. 2022).We align the images to Gaia DR3 (Gaia Collaboration et al. 2021), co-add, and finally drizzle them (Fruchter & Hook 2002) to a 0. 02 pixel scale for the Short Wavelength (SW) NIRCam bands and to 0. 04 1 Supplementary material and catalogs of selected sources: 10.5281/zenodo.7614908;mosaics and field catalogs 10.17894/ucph.e3d897af-233a-4f01-a893-7b0fad1f66c2 2 10.5281/zenodo.7143382for all the remaining JWST and Hubble Space Telescope (HST ) filters.We provide further details about the individual fields in Appendix A.

Extraction
We extract sources using a detection image produced by combining of all the "wide" (W) NIR-Cam Long Wavelength (LW) filters available (typically F277W+F356W+F444W) optimally weighted by their noise maps.For source extraction, we use sep (Barbary 2016), a pythonic version of source extractor (Bertin & Arnouts 1996).We extract the photometry in circular apertures with a diameter of 0. 5 and correct to the "total" values within an elliptical Kron aperture (Kron 1980) 3 .The aperture correction is computed on the LW detection image and applied to all bands.The depth in the reference 0. 5 apertures in the five NIRCam bands that we require to select candidate quiescent galaxies (F150W, F200W, F277W, F356W, and F444W, Section 3) are reported in Table 1.The galaxy distribution as a function of redshift for F444W is shown in Appendix (Figure A.1) and for the remaining bands in Figure Set A1.An extra-correction of ∼ 10% to account for the flux outside the Kron aperture in HST bands and optimal for point-like sources is computed by analyzing the curve of growth of point spread functions (PSF)4 .

Modeling of the spectral energy distribution
We utilize eazy-py5 (Brammer et al. 2008) to estimate photometric redshifts, rest-frame colors, and stellar masses from the 0. 5 diameter aperture photometry corrected to total fluxes as described above.We apply residual zeropoint corrections to optimize the photometric redshifts with solutions free to vary in the interval z = 0 − 18.We use the same set of 13 templates from the Flexible Stellar Populations Synthesis code (fsps, Conroy & Gunn 2010) described in Kokorev et al. (2022) and Gould et al. (2023), linearly combined to allow for the maximum flexibility.This set of templates covers a large interval in ages, dust attenuation, and log-normal star formation histories -spanning the whole U V J rest-frame color diagram.More specifically, the corr sfhz 13 subset of models within eazy contains redshift-dependent SFHs, which, at a given redshift, ex-Figure 1. U − V , V − J rest-frame color diagrams for our combined sample of galaxies in JWST fields binned in redshift as labeled.Filled circles indicate our visually-inspected U V J quiescent sample and their 1σ uncertainties, color coded according to their stellar mass.We circled in black the object in the "strict" sample.Gray points indicate the rest of the sample at those redshifts (Section 3).The color intensity scales as the density of points.The red dotted and solid lines indicate the standard selection box (Williams et al. 2009) and a looser version allowing for an extra pad of 0.2 mag, respectively.The black arrow shows the effect of reddening for AV = 1.clude histories that start earlier than the age of the Universe.A template derived from the NIRSpec spectrum of a confirmed strong line emitter at z = 8.5 -ID4590 from Carnall et al. (2022a) -is also included to allow for an extra degree of freedom in photometric solutions of distant objects, but not accounted for in the stellar mass calculation6 .The templates are created adopting a Chabrier (2003) initial mass function and applying a Kriek & Conroy (2013) dust attenuation law (dust index δ = −0.1,R V = 3.1), where the maximum allowed attenuation is also redshift-dependent.A fixed grid of nebular emission lines and continuum from cloudy (v13.03)models is added to the templates within fsps (metallicity: log(Z/Z ) ∈ [−1.2, 0], ionization parameter log(U ) = −1.64,−2; Byler 2018).Given their fixed ratios and sole purpose of modeling the photometry, the strength of the emission lines should be taken with caution.We, thus, do not include them in our analysis.The templates, their input parameters, and the redshift evolution of their allowed SFHs and attenuation are available online 7 .Also, a correction for the effect of dust in the Milky Way is applied to the templates within eazy-py, pulling the Galactic dust map by Schlafly & Finkbeiner (2011) from dustmaps (Green 2018).This effect is relevant for SMACS0723 (E(B − V ) MW = 0.19, see also Faisst et al. 2022) and to a lesser extent for the rest of the fields (E(B −V ) MW = 0.007−0.07).In terms https://github.com/gbrammer/eazy-photoz/tree/master/templates/sfhz of photometric redshifts, we obtain a good agreement with spectroscopic determinations from archival observations when available (G. Brammer et al., in preparation).For reference, in the CEERS field we estimate a σ NMAD = 0.0268 (0.0187) for the spectroscopic sample (clipping catastrophic outliers), respectively8 .Stellar masses are consistent with independent estimates obtained with finer grid-based codes (Figure B.3 in Appendix for a comparison with 3D-HST ).Star formation rates are also found in agreement with determinations with alternative codes at z = 0.5−3 (Gould et al. 2023).However, we opt not to rely on SFR for our selection and analysis at z > 3 to adhere as close as possible to observables.

Rest-frame colors
As described in Gould et al. (2023), besides photometric redshifts constrained by spectral features, eazy-py returns the physical quantities attached to each template, propagated through the minimization process and computed using the same set of coefficients providing the best-fit z phot .Uncertainties are estimated at z phot as the 16-84% percentiles of a 100 fits drawn from the best-fit template error function.We compute the restframe magnitudes in the GALEX N U V band (λ = 2800 Å), the U and V Johnson filters defined as in Maíz Apellániz (2006), and the J 2MASS passband (Skrutskie et al. 2006).The rest-frame magnitudes are computed following a hybrid approach that uses the templates as guides for a weighted interpolation of the observations and accounts for bandpasses and relative depths (Brammer et al. 2008, 2011, and Appendix A in Gould et al. 2023).This allows for a color determination which relieves the dependency on the adopted models, while using the whole photometric information.

SAMPLE SELECTION
Before selecting quiescent galaxy candidates at z > 3, we start by applying a series of loose cuts to immediately reject galaxies with unreliable photometric modeling.We constrain the quality of the fit to χ 2 /N filt ≤ 8, where N filt ≥ 6 is the number of available filters.The latter includes NIRCam wide bands at 1.5, 2.0, 2.7, 3.5, 4.4 µm in every field with the exception of SPT2147, where imaging with F150W was not taken.Coupled with the adoption of the NIRCam LW detection image (Section 2.1), this requirement enforces a selection based on JWST data, while the coverage at observed wavelengths longer than 3 µm allows for robust determinations of stellar masses.Then we apply a loose cut at > 5 on the quadratic sum of the S/N of the aperture fluxes in these NIRCam bands.We constrain the location of the peak and the tightness of the redshift probability distribution function p(z) (max{p(z)} > 0.5, (p(z) 84% − p(z) 16% )/(2 p(z) 50% ) < 0.3, where p(z) i% indicate the i-th percentile of p(z)).Finally, we apply a cut in redshift at 3 ≤ z ≤ 6.5 with a buffer of dz = 0.1 and accounting for the uncertainty on the best-fit solution (p(z) 84% ≥ 3 − dz and p(z) 16% ≤ 6.5 + dz).To pick quiescent objects, we opt for a rest-frame color selection following a dual approach.

U V J color diagram
On the one hand, we select objects in the classical U V J diagram.This allows us to directly compare our results with a large body of literature that has accumulated over the last two decades.We allow for a 0.2 mag extra pad when compared with the cuts in Williams et al. (2009) and we initially retain sources with 1σ uncertainties on the colors consistent with the selection box as long as σ color < 0.5 mag.We then visually inspect the images and the SED fits of 251 candidate quiescent galaxies.We retain 109/251 objects (∼ 45%) after excluding remaining bad fits or poor quality images affected by edge effects, spikes, or contaminating bright sources.We show the location of the visually inspected sample in three redshift bins at z > 3 in Figure 1.The visual selection significantly shrinks the initial pool of candidates.This is expected given the deliberate choice of starting from rather loose constraints not to lose potential good candidates.The visual cut particularly hits the highest redshift pool of candidates: we retain 3/56 galaxies at z > 5 largely due to poor quality SEDs.For transparency, all of the SEDs and cutouts of the discarded sources during the visual check are also released.
To draw straightforward comparisons with previous works and in attempt to remedy the larger contamination that inevitably affects our expanded selection box, we further flag our sources as "strict" and "padded".The first tag refers to 55/109 sources that fall in the classical QG box, also accounting for their 1σ uncertainties (34/109 without including the latter as in the "standard" selection).The second flag refers to 67/109 sources that have nominal (i.e., without including uncertainties) colors within the 0.2 mag padded locus of QGs.The overlapping sample comprises 51 galaxies.Differences in the derived number densities primarily reflect these further distinctions.
Three-color NIRCam SW and LW cutouts, photometry, SED models, and basic properties estimated with eazy-py are publicly available for the U V J-selected samples.
In parallel, we follow the novel method described in Gould et al. (2023) (see also Antwi-Danso et al. 2022 for an alternative approach introducing a synthetic band).The authors incorporate the N U V magnitude in their selection and model the galaxy distribution in the (N U V − U , U − V , V − J) space with a minimal number of Gaussians carrying information (Gaussian Mixture Model, GMM, Pedregosa et al. 2011).The addition of the N U V magnitude makes the selection more sensitive to recent star formation and, thus, to recently quenched or post-starburst objects (Arnouts et al. 2013;Leja et al. 2019b), which are expected to be observed at high redshift as we approach the epoch of quenching of the first galaxies (D'Eugenio et al. 2020;Forrest et al. 2020b).Moreover, the GMM allows to fully account for the blurred separation between starforming and quiescent galaxies at z > 3, assigning a "probability of being quiescent" P Q to each object and bypassing the use of arbitrary color cuts.The GMM grid is calibrated on a sample of 2 < z < 3 galaxies in the COSMOS2020 catalog (Weaver et al. 2022a) assuming 5× more conservative Spitzer /IRAC uncertainties and refit with eazy-py in a similar configuration to that adopted here (Valentino et al. 2022).To account for the uncertainties on the colors, we bootstrap their values 1000 times and use the median and 16-84% percentiles of the distribution as our reference P Q,50%9 and its 1σ uncertainties.We also list the nominal P Q associated with the best-fit colors in our catalogs.
In the rest of our analysis, we adopt a cut at P Q,50% ≥ 0.1 to select candidate quiescent galaxies, with a threshold set at P Q,50% = 0.7 to separate passive galaxies from objects showing features compatible with more recent quenching (see Gould et al. 2023 for a description of the performances of different cuts benchmarked against simulations).As for the U V J-selected sample, we visually inspect all of the images and SEDs of the candidates that made our initial P Q,50% cut.Finally, we retain 50/71 sources (70%) with 0.1 ≤ P Q,50% < 0.7 and 18/20 (∼ 90%) truly passive galaxy candidates with P Q,50% ≥ 0.7.Their location in the projected N U V −U , V − J plane is shown in Figure 2.

Overlap between selections
As noted in Gould et al. (2023), a selection in the N U V U V J arguably outperforms the classical U V J in selecting quiescent (passive and recently quenched or post-starbust) galaxies at z > 3.However, the two criteria partially overlap and identify the same quiescent sources -to an extent fixed by the P Q threshold and the exact location of the selection box in the U V J diagram.The boundaries adopted here slightly differ from those in Gould et al. (2023), but the resulting overlap between the selection criteria at 3 < z < 5 is similar.Focusing on the visually inspected samples, 52 sources are selected by both techniques.These amount to ∼ 75% and ∼ 50% of the N U V U V J and U V Jselected objects, respectively, comparable with the fractions reported in Gould et al. (2023) in the same redshift range.In more detail, 100% and ∼ 70% of the sources with P Q,50% ≥ 0.7 and 0.1 ≤ P Q,50% < 0.7 are part of the U V J sample.Moreover, 16/18 objects with P Q,50% ≥ 0.7 fall within the standard U V J selection box.This is because the galaxies assigned lower P Q,50% values are those in the region bordering star-forming and quiescent, whereas galaxies with higher P Q,50% are those which resemble more classically quiescent galaxies owing to how the model is trained (Figure 2).Therefore, it is expected that the U V J selected sample has a smaller overlap with galaxies that have lower P Q,50% values and the galaxies that it does not select are those which are recently quenched.The overlap is reflected also in the M distributions of the selected samples (Figure C.7 in Appendix).Lower P Q,50% values are associated to bluer, more recently quenched, but also lower mass objects, otherwise missed by U V J selections.

Sanity checks on the sample
We test our selection and draw comparisons with what has been achieved before the advent of JWST in a variety of ways, as summarized below.More details can be found in Appendix B.

Comparison with HST-based photometry
First, we compare our HST /F160W photometry (consistent with that of NIRCam/F150W), photometric redshift, and stellar mass estimates against those from the 3D-HST catalog (Skelton et al. 2014) overlapping with part of the CEERS (EGS) and PRIMER (UDS) areas (Figure B.3).Despite different detection images, for those sources in common we find excellent agreement in terms of aperture photometry and redshifts derived.Moreover, minimal systematic offsets in log(M /M ) (< 0.2 dex) make our results fully consistent between different SED modeling codes (Appendix B.1).

Availability of HST imaging
We also test the impact on our sample selection of HST filters, which increase the sampling of the restframe UV/optical wavelengths in some of the fields (Table 3 in Appendix A).We refit the photometry in the CEERS and PRIMER fields retaining only the available NIRCam filters among those at 0.9, 1.15, 1.5, 2.0, 2.7, 3.5, and 4.4 µm, mimicking the situation for Stephan's Quintet where no HST imaging is at disposal.We obtain fully consistent z phot , M , and rest-frame color estimates within the uncertainties, especially when removing the effect of z phot from the calculation and focusing on the 3 ≤ z ≤ 6.5 interval of interest (Figure B.4).This holds also when F090W, probing wavelengths shorter than N U V at the lower end of the redshift interval that we explore, is not available, as in the case of CEERS.We also re-apply the N U V U V J selection, including bootstrapping, and obtain consistent samples taking into account the uncertainties.

Dusty star-forming or high-redshift contaminants
When available, we look for counterparts at long wavelengths to exclude obvious dusty interlopers.We search for matches in sub-millimetric observations in CEERS (450 and 850 µm with Scuba-2, Zavala et al. 2017;Geach et al. 2017), PRIMER (870 µm with ALMA from the AS2UDS survey, Dudzevičiūtė et al. 2020, andCheng et al. 2022), and SMACS0723 (1.1 mm with ALMA from the ALCS Survey, Kokorev et al. 2022;S. Fujimoto et al. in preparation).We found only one potential association with a ∼ 5σ Scuba-2 detection at 850 µm in CEERS (S2CLS-EGS-850.063 in Zavala et al. 2017, #9329 in our catalog).This is a weak sub-millimeter galaxy candidate possibly associated with an overdensity (S.Gillman et al. in preparation).Removing it from our sample does not change the results of this work.The absence of longwavelength emission in our pool of visually inspected candidates supports their robustness.We also note that the contamination of high-redshift Lyman Break Galaxies is negligible, given their number densities (Fujimoto et al. 2022).

JWST-based photometric samples
We correctly identify and select the spectroscopically confirmed QGs in Schreiber et al. (2018b) at consistent z phot .Our relaxed U V J cuts and P Q,50% ≥ 0.1 recover 14/17 and 11/17 candidate QGs selected via SED modeling and a sSFR threshold in Carnall et al. (2022b), respectively.When considering only their "robust" sample, there is a 89% and 78% overlap with our U V J and P Q,50% selections.We calculate lower, but generally consistent z phot (Figure B.5, see also Kocevski et al. 2022).Minor systematic offsets in the M and F200W magnitudes are present (Figure B.5).

NUMBER DENSITIES
For each field, we compute the comoving number density n of candidate quiescent galaxies in three redshift bins (z ∈ [3, 4), [4,5), and [5, 6.5)).We compute the number of objects per bin by integrating the p(z), thus accounting for the uncertainties on the photometric redshift determination.As an alternative estimate of the statistical errors associated with the latter, we randomly extract 1000× each z phot within their p(z) and compute the median and 16 − 84% percentiles, finding consistent results.We also compute the 68% Poissonian confidence intervals or upper limits (Gehrels 1986).The comoving volumes are calculated starting from the area subtended by the observations and satisfying the requirements in terms of band coverage and minimum number of filters (Section 3).For the Sunrise and SMACS0723 fields, we account for the effect of gravitational lensing on the volume at z = 3 − 5 as in Fujimoto et al. (2016).We estimate the intrinsic survey volume by producing magnification maps at z = 3, 4, and 5.We base the calculation on the mass model constructed with the updated version of glafic (Oguri 2010(Oguri , 2021) ) using the available HST and JWST data (Harikane et al. 2022;Welch et al. 2022a).The effect of lensing on the effective area varies negligibly in the redshift interval and luminosity regime spanned by our samples of candidate quiescent galaxies.The linear magnification for the only QG candidate in proximity of SMACS0723 is ∼ 1.7.For two candidates in WHL0137 (Sunrise) with log(M /M ) < 9.5 and P Q,50% < 0.7, this factor is ∼ 3. We did not apply the magnification correction to the parameter estimates.This does not affect the conclusions on number densities.
In Figures 3 and 4, we show the p(z) and the corresponding comoving number densities for the U V J and N U V U V J-selected quiescent galaxies in each of the 11 fields that we consider.A combined estimate based on the aggregate area of 145.1 arcmin 2 is also presented.In each panel, we report the number densities in stellar mass bins of 10 9.5 ≤ M < 10 10.6 M and M ≥ 10 10.6 M .The high mass threshold is chosen to directly compare these results with those in the litera-Figure 3. Number densities of U V J-selected galaxies in public JWST fields.The purple and green colors mark the p(z) of individual robust "strict" U V J quiescent candidates with M ≥ 10 10.6 M and 10 9.5 ≤ M < 10 10.6 M , respectively.The p(z) are normalized by their area ( z p(z) dz = 1).The sky coverage of each field is as labeled.We included the effect of gravitational lensing at these redshifts in the calculation of the areas around galaxy clusters (Sunrise, SMACS0723), for which the masses should be intended as observed (not delensed).The comoving number densities in units of 10 −5 Mpc −3 obtained from the integration of p(z) within the 3 < z < 4, 4 < z < 5, and 5 < z < 6.5 bins marked by dotted lines are reported.The errors mark the 68% Poissonian confidence intervals.Estimates of the uncertainties from bootstrapping are within brackets.The first and second rows indicate n in the 10 9.5 ≤ M < 10 10.6 M and M ≥ 10 10.6 M bins, respectively.
3 < z < 4 [9.5, 10.6) 3.9 +1.2 −0.9 4.1 +1.2 −0.9 Upper limits are at 1σ using the same approach (Gehrels 1986).Statistical uncertainties are accounted by integrating the p(z) within the redshift intervals.The uncertainties due to cosmic variance are expressed as fractional σCV deviations (Section 4.1).The selections are described in Section 3. The adopted threshold for the ture (Table 4).Such a threshold also allows us to safely compare different fields.This is clear from Figure C.6 in Appendix C, showing the stellar mass limit in each field for the overall sample of galaxies and QGs.The number density estimates that we derive for the combined field are reported in Table 2.
For the U V J-selected galaxies, we show the results for the sources "strictly" obeying the classical selection, while accounting for the color uncertainties.Adopting the "padded" sample returns consistent results, while estimates based on the whole "robust" pool of galaxies would be intended as an upper limit.An even stricter criterion accounting only for galaxies with nominal color within the standard U V J color box returns 2× and 1.2× lower, but fully consistent number densities in the lower and higher mass bins, respectively.In principle, an Eddington-like bias could be introduced by our "strict" selection coupled with the asymmetric distribution of galaxies in the color and mass space (more blue star-forming and lower mass systems can scatter into the selection box than red massive quiescent candidates that move out).However, this effect seems to be of the same order of magnitude of the statistical uncertainties.We also note that we conform to a pure color selection and we do not apply any formal correction for contamination of dusty interlopers (∼ 20% for standard U V J selection, Schreiber et al. 2018b, thus likely higher for the padded and the robust samples).The N U V U V Jselected sample (P Q,50% ≥ 0.1) is as numerous as the U V J one, despite the partial overlap between the two criteria (Section 3.3), similar to what is found in Gould et al. (2023).

Cosmic variance
To compute the cosmic variance, we use the prescription of Steinhardt et al. (2021), which is based on the cookbook by Moster et al. (2011).These authors assume a single bias parameter that links stellar to halo masses in log(M /M ) bins of 0.5 dex.In simulations, this assumption has been found to be valid only to a 0.2 dex level even for massive galaxies (Jespersen et al. 2022;Chuang & Lin 2022).However, this is < 1/5 of the bin widths used in this paper and, thus, the approximation should be appropriate.The cosmic variance is computed for each individual field taking into account the survey geometry.In order to get the cosmic variance for the bin at log(M /M ) = 9.5 − 10.6, we weight the contribution of each 0.5 dex bin to the total cosmic variance by the relative number densities in Weaver et al. (2022b).The total cosmic variance is then computed as: which assumes that all fields are independent.10 .The relative uncertainties due to cosmic variance in the combined field are reported in Table 2.

4.2.
A compendium of number densities of massive quiescent galaxies at 3 < z < 4 Excellent depictions of how many QGs each individual survey or simulation find as a function of redshift are available in the literature (e.g., Straatman et al. 2014;Merlin et al. 2019;Girelli et al. 2019;Cecchi et al. 2019;Shahidi et al. 2020;Weaver et al. 2022b;Gould et al. 2023;Casey et al. 2022;Long et al. 2022).However, drawing direct comparisons among different works and evaluating the impact of various selections is complicated by the introduction of systematic assumptions.In Figure 5, we attempt to partially remedy this situation by reporting number densities at least adopting a consistent redshift interval (3 < z < 4) and lower mass limit for the integration (10 10.6 M ) for similar IMFs (Chabrier 2003;Kroupa 2001).We also add number density estimates from EAGLE (Schaye et al. 2015;Crain et al. 2015), Illustris (Vogelsberger et al. 2014), Illustris-TNG 100, and 300 simulations (Nelson et al. 2018).We count simulated galaxies with sSFR ≤ 10 −10 yr −1 within 4× and 2× the half-mass radius for EA-GLE and Illustris(-TNG), respectively (see Donnari et al. 2019 for a discussion on different selection criteria of QGs, average timescales to estimate SFRs, and physical apertures in simulations).We consider snapshots at z = 3.0 and z = 3.7−3.9(Valentino et al. 2020).
Our number density estimates from the combined fields are of the order of ∼ 2.5 × 10 −5 Mpc −3 , consistent with some of the determinations with similar color or sSFR cuts (Schreiber et al. 2018b;Merlin et al. 2019).Our estimates are ∼ 2× larger than the most recent measurements in the largest contiguous survey among those considered, COSMOS (Weaver et al. 2022b), also when adopting very consistent color selections (Gould et al. 2023).Interestingly, earlier determinations in the same field retrieved significantly lower estimates (Muzzin et al. 2013;Davidzon et al. 2017;Girelli et al. 2019;Cecchi et al. 2019).This is due to a combination of deeper and homogeneous measurements in the optical and near-IR over a twice larger effective area, now available in COSMOS2020 (Weaver et al. 2022a), more conservative and pure samples of QGs, the specific templates used in each work, and the integration of best-fit Schechter function underestimating the observed values at the high end of the stellar mass function.
The new detection and extraction based on JWST LW observations allows for the selection of redder sources and improved deblending that was previously based on HST bands or ground-based observations.This allows for a better identification of higher-redshift, less massive quiescent galaxies, and more robust M by finding breaks at longer wavelengths, pinpointing objects with lower mass-to-light ratios, and removing blended objects (see also the discussion in Carnall et al. 2022b).Nevertheless, in our most massive bin, the variations among different works are still dominated by systematics in the selection, modeling, and cosmic variance.

Field variations and groups of quiescent galaxies
We notice a substantial field-to-field variation especially when focusing on the most massive galaxies.Compared with the average number density on the full combined field of ∼ 145 arcmin 2 , we find per field value oscillations of a factor of 2 − 3× (Table 5).We ascribe these differences to cosmic variance and to the fact that massive quiescent systems might already be signpost of distant overdensities and protoclusters − massive halos able to fast-track galaxy evolution.In Table 5, we report the fractional 1σ uncertainties due to cosmic variance in each field.Taken individually, an uncertainty of ∼ 30 − 50% affects the number densities in the two mass bins for the largest contiguous areas that we considered.The Stephan's Quintet and CEERS fields are emblematic in this sense, appearing under-and over-dense despite a similar sky coverage.CEERS displays the largest number densities of QG > 10 10.6 M in our compilation, consistent with the estimates in Carnall et al. (2022b).There we find a remarkable pair of candidate quiescent systems with consistent z phot = 3.54−3.38(#9622-9621, also "robust" and not "robust" candidates in Carnall et al. 2022b, respectively).The pair, possibly interacting, is surrounded by two more red sources with similar z phot ∼ 3.18 − 3.54 that fall in the visually vetted U V J sample (#9490, 9329).This is reminiscent of the massive galaxies populating the "red sequence" in clusters, also used to find evolved protostructures at high redshift (Strazzullo et al. 2015;Ito et al. 2023).Similar QGs have already been found in overdensities at z 3 (Mc-Conachie et al. 2022;Kalita et al. 2021;Kubo et al. 2021) or in close pairs with other massive objects ("Jekyll & Hyde" at z = 3.717, Glazebrook et al. 2017;Schreiber et al. 2018a, of which a pair of quiescent objects would be a natural descendant).Two of these pairs or small collections of red galaxies at similar redshifts and close in projection are in our list -not surprisingly, especially in the loosest U V J-selected sample.

A look at lower masses
Figures 3 and 4 show the number densities at lower stellar masses (10 9.5 ≤ M < 10 10.6 M ), now safely accessible with JWST also at such high redshifts.In fact, the low-mass end of the M distributions starts declining at thresholds as low as ∼ 10 8 M at 3 < z < 6.5 (excluding SPT2147, Figure C.6 in Appendix C).The lower limit of M = 10 9.5 M chosen for the calculation is similar to that fixed in some of the works listed in Table 4. Thus, it allows us to derive a relatively straightforward comparison, discounting some of the systematics mentioned above.At 3 < z < 4, we estimate modestly higher (∼ 1.2 − 1.6×) number densities than in the most massive bin, but in agreement within the uncertainties.The difference with previous works is similar in each bin, when available.This is consistent with the expected shape of the stellar mass function of red quiescent galaxies, roughly peaking and flattening or turning over at ∼ 10 10.6 M at these redshifts (Weaver et al. 2022b) and revealing a steady build up of lower-mass QGs (Santini et al. 2022).Promising low-mass QGs have been already confirmed with JWST /NIRISS at z ∼ 2.5 in the GLASS field (Marchesini et al. 2023).We defer to future work a comprehensive analysis of the stellar mass functions of quiescent populations at these redshifts.

High-redshift candidates
We can now look at galaxies at 4 < z < 5.As in the 3 < z < 4 interval, when considering the combined fields, we estimate number densities that are ∼ 2.5−4.5×smaller in above and below M = 10 10.6 M than those in CEERS by Carnall et al. (2022b), who also perform a selection starting from JWST images.The difference shrinks to a factor of ∼ 2.5−1× when we consider only the same field and the "robust" sample in their work.This seems to suggest that cosmic variance and early overdense environments are effective in producing substantial field-to-field variations also at z > 4 (Section 4.3).For reference, our number density estimates in the same massive bin (M ≥ 10 10.6 M ) are consistent with those in the COSMOS field by Weaver et al. (2022b), but 1.8× larger than what retrieved in the same field but using a color selection similar to ours (Gould et al. 2023; see the discussion therein on the agreement with the latest COSMOS2020 number densities).When integrating down to lower mass limits (10 9.5 M , but not homogenized among different works at this stage, given the impact of different depths at z > 4), we retrieve similar n as in COSMOS ((1.0 ± 0.3) × 10 −5 Mpc −3 for M > 10 9.9 M , Weaver et al. 2022b) and 1.6 − 2× larger than in large-field HST surveys such as CAN-DELS (∼ 7.9 × 10 −6 Mpc −3 for the "complete" sample at M > 5 × 10 9 M , Merlin et al. 2019).Finally, the upper limits on number densities for the highest redshift bin at 5 < z < 6.5 should be taken with caution, given the area covered in our analysis (Table 2).
Focusing on the highest envelope of the redshift interval spanned by our sample, we find a few promising candidates at z > 4.5.The SEDs and three-color images of the 5 most robust sources falling either in the "strict" or "padded" U V J selections are shown in Figure 6.We do not find reliable candidates at z 5.2, which signposts the earliest epoch of appearance of quiescent objects in our current sample 11 .Source #185 in PRIMER (z = 4.51 +0.16  −0.26 , log(M /M ) = 10.9) is also picked as among the most reliable quiescent candidates by the N U V U V J criterion (P Q,50% = 0.79).An entry with z ∼ 3.2 at < 0. 3 from this source is present in previous catalogs of this field (Skelton et al. 2014;Mehta et al. 2018), but more consistent and blended with the nearby blue object (a chance projection in our analysis, z phot = 2.78 +0.10 −0.08 ).The remaining sources are assigned lower P Q,50% values compatibly with their bluer colors and more recent or possibly ongoing quenching.All these candidates appear rather compact.Sources # 789 and 303 in PRIMER are compatible with the locus of stars in the FLUX RADIUS, MAGAUTO plane and should be taken with a grain of salt.For comparison, we also checked the "robust" objects at z > 4.5 in Carnall et al. (2022b).However, we retrieve only # 101962 (our #2876) above z = 4 (Appendix B.4, see also Kocevski et al. 2022).Direct spectroscopic observations with JWST are necessary to break the current ceiling at z ∼ 4 imposed by atmospheric hindering to groundbased telescopes and confirm the exact redshifts of these candidates.

Revisiting the comparison with simulations
11 Two bluer and potentially quenching objects are picked by the loosest U V J selection at z > 5. Their SEDs and color images are part of the overall release.
For what concerns simulations, if we limit our conclusions to the homogenized massive bin and 3 < z < 4 interval, we find a broad agreement with the Illustris TNG suite at the lower end of the redshift range (z = 3) and a rapidly increasing tension above this threshold (Valentino et al. 2020).EAGLE and Illustris seem to struggle to produce massive M ≥ 10 10.6 M QGs already at z = 3, while the situation seems partially alleviated if one includes lower mass galaxies in the calculation (Merlin et al. 2019;Lovell et al. 2022).Spectroscopically confirmed massive QG at z > 4 would not only exacerbate the tension with these simulations, but also for the latest-generation examples, such as FLARES (Lovell et al. 2021;Vijayan et al. 2021).We do not find any M > 10 10.6 M objects with sSFR < 10 −10 yr −1 in EAGLE or Illustris-TNG at z = 4 and above, while FLARES produces 2 dex lower number densities at z = 5 (n = 7.2 × 10 −8 Mpc −3 , Lovell et al. 2022).

CONCLUSIONS
We present a sample of ∼ 80 JWST -selected candidate quiescent and quenching galaxies at z > 3 in 11 separate fields with publicly available imaging collected during the first 3 months of telescope operations.We homogeneously reduce the JWST data and combine them with available HST optical observations.We both perform a classical U V J selection and apply a novel technique based on Gaussian modeling of multiple colors -including an N U V band sensitive to recent star formation, which is necessary to explore the quenching of galaxies in the early Universe.Here we focus on a basic test for simulations and empirical models: the estimate of comoving number densities of this population.
• We estimate n ∼ 2.5 × 10 −5 Mpc −3 for massive candidates (≥ 10 10.6 M ) with both selections, but substantial field-to-field variations of the order of 2 − 3×.This is likely due to cosmic variance (∼ 30 − 50% uncertainty in the largest contiguous fields of ∼ 30 arcmin 2 such as CEERS or Stephan's Quintet) and the fact that early and evolved galaxies might well trace matter overdensities and the emerging cores of protoclusters already at z > 3. We find promising candidate pairs or groups of quiescent or quenching galaxies with consistent redshifts in the field with the highest number density.
• We compile and homogenize the results of similar attempts to quantify the number densities of massive QG at 3 < z < 4 in the literature.The comparison across almost 20× different determinations highlights the impact of cosmic variance and systematics primarily in the selection techniques.
The most recent estimates seem to converge toward a value of n ∼ 1 − 2 × 10 −5 Mpc −3 -not exceedingly far from what established via groundbased observations.
• We apply our homogenization also to publicly accessible large-box cosmological simulations.As noted in the literature, a tension with observations at increasing redshifts is evident -to the point that even a single confirmation of a massive QG at z > 4 − 4.5 would challenge some of the theoretical predictions.A handful of promising candidates up to z ∼ 5 are found in our systematic search and presented here.
• We start exploring the realm of lower mass QG candidates, taking advantage of the depth and resolution of JWST at near-IR wavelengths.We measure number densities at 10 9.5 ≤ M < 10 10.6 M similar to those at ≥ 10 10.6 M , consistent with the expected flattening or turnover of the stellar mass function of quiescent objects and the onset of the low-mass quenched population.
This work is the first of a series of articles that will focus on the characterization of several aspects of the sample selected here (morphologies, SED and SFH modeling, also resolved, and environment).We remark that all of the high-level science products (notably We acknowledge the careful reading and the constructive comments from the anonymous referee.We warmly thank Emiliano Merlin, Giacomo Girelli, Abtin Shahidi, and Micol Bolzonella for computing and sharing their number densities in the specified redshift and mass intervals.We also thank Dan Coe for sharing the magnification maps computed by the RELICS team.This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope.The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST.The specific observations analyzed can be accessed via 10.17909/g3nt-a370.These observations are associated with programs ERS #1324, 1345, and 1355; ERO #2736; GO #1837 and 2822; GTO #2738; and COM #1063.The authors acknowledge the teams and PIs for developing their observing program with a zero-exclusive-access period.The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140.S.F. acknowledges the support from NASA through the NASA Hubble Fellowship grant HST-HF2-51505.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.M.H. acknowledges funding from the Swiss National Science Foundation (SNF) via a PRIMA Grant PR00P2 193577 "From cosmic dawn to high noon: the role of black holes for young galaxies".This work was supported by JSPS KAKENHI Grant Numbers JP21K03622, 20K14530, and 21H044902.K.I. acknowledges support from JSPS grant 22J00495.G.E.M. acknowledges the Villum Fonden research grants 13160 and 37440.O.I. acknowledges the funding of the French Agence Nationale de la Recherche for the project iM-AGE (grant ANR-22-CE31-0007).(Fruchter & Hook 2002;Koekemoer et al. 2003), glafic (Oguri 2010(Oguri , 2021)).

APPENDIX
A. FIELDS Here we provide a brief summary of the available observations for each field.The JWST and HST imaging availability is summarized in Table 3 (Bagley et al. 2022 for a full description of the program and an official data release of the CEERS team).In this work, we made use of the NIRCam imaging in the "wide" F115W, F150W, F200W, F277W, F356W, and F444W filters, plus the "medium" F410M band.We incorporated available HST observations from the archive (CHArGE, Kokorev et al. 2022).

A.2. Stephan's Quintet
Stephan's Quintet has been targeted and the images immediately released as part of the Early Release Observations (ERO # 2736, Pontoppidan et al. 2022).No HST coverage is available in our archive.In Figure A.2, we show the nominal overlap of the filters that we required for the selection, but we carved a large portion of the central part of the field where contamination from the galaxies belonging to the group was too high to ensure good quality photometry.This effectively reduced the area by ∼ 5 arcmin 2 .

A.3. PRIMER
The Public Release IMaging for Extragalactic Research (PRIMER, GO 1837, PI: J. Dunlop) is a Cycle 1 accepted program targeting contiguous areas in the COSMOS (Scoville et al. 2007) and Ultra-Deep Survey (UDS, Lawrence et al. 2007) fields with NIRCam and  Treu, Treu et al. 2022).The parallel fields are sufficiently far from the cluster that gravitational lensing does not appreciably affect our work.Abell 2744 has been targeted by several HST programs, including the Grism Lens-Amplified Survey from Space (GLASS) itself and a project tailored to maximally exploit the scientific return of the parallel fields (GO/DD 17231, PI: T. Treu), which we also included in our data.

A.7. Sunrise
We dubbed the cluster-lensed field WHL0137-08 from the Reionization Lensing Cluster Survey (RELICS, Coe et al. 2019, GO #2822) after the "Sunrise arc" that was discovered in it (Salmon et al. 2020;Vanzella et al. 2022) -even hosting a highly magnified star (Welch et al. 2022b).Being part of RELICS, ample HST ancillary data is available.The area reported in Table 1 accounts for the lensing effect at z = 3 − 5. JWST data were included in updated magnification maps generated with glafic (Oguri 2010(Oguri , 2021)).

A.8. SMACS0723
SMACS0723 is also part of the RELICS survey, one of the spectacular, and immediately released ERO objects (Pontoppidan et al. 2022).Here we made use of NIRCam imaging from detectors targeting the cluster and a position offset from it.MIRI, when available, was included.Also in this case we accounted for the effect of lensing in the area centered on the cluster using an updated version of previous magnification maps with glafic, now including JWST data.HST coverage from RELICS is available.Long-wavelength observations from the ALMA Lensing Cluster Survey (PI: K. Kohno) were used to look for possible dusty contaminants, when available.
A.9. SGAS1723, SPT0418, and SPT2147 These are fields from the Targeting Extremely Magnified Panchromatic Lensed Arcs and Their Extended Star formation (TEMPLATES) ERS program (ERS 1355, PI: J. Rigby).The primary targets are 4 strongly galaxy-lensed systems with ample ancillary data across the electromagnetic spectrum.In this work, we relied on the imaging portion of the ERS program.Single galaxy lensing does not affect the field on large scales.SPT2147 was not imaged with the F115W and F150W filters.

B. SANITY CHECKS ON THE SAMPLE
Here we describe in more detail the tests on the robustness of our sample selection to which we briefly referred in Section 3.4.

B.1. Photometry and source extraction
We compared our photometric extraction and SED modeling (Section 2) with those from 3D-HST (Skelton et al. 2014) for sources in CEERS (EGS) and PRIMER (UDS).We matched sources allowing for a maximal < 0. 5 separation.Figure B.3 shows the comparison in z phot , M , total HST /F160W photometry computed from our reference 0. 5 diameter aperture, and that within a common 0. 7 aperture.The agreement is overall excellent, despite different detection images and corrections applied.The aperture photometry is fully consistent and so are the z phot estimates from the previous and current version of eazy(-py).The F160W total magnitudes computed starting from the 0. 5 diameter apertures considered here are fainter than those from 0. 7 in 3D-HST in CEERS and PRIMER: the median differences are 0.11 (σ MAD = 0.14) and 0.16 (0.14) mag, respectively.However, at fixed aperture, the total magnitudes are fully consistent.The difference arises from the detection bands and where the aperture correction is computed: F160W for 3D-HST and the NIRCam combined LW image in our analysis.Finally, the total M are systematically lower in 3D-HST than in our JWST -based catalogs of CEERS and PRIMER, with median differences and MAD of 0.19 (σ MAD = 0.30) and 0.13 (0.24) dex, respectively.All things considered, the offsets are fully ascribable to the different recipes adopted to estimate these quantities and consistent with typical systematic uncertainties inevitably present when we compare different catalogs of the same sources.Our samples of U V J and N U V U V J selected QG at z > 3 do not appreciably deviate from these trends.

B.2. Availability of HST photometry
We tested our sample selection against the availability of HST filters sampling the rest-frame UV/optical emission at z > 3.As mentioned in Section 3.4, we refitted the photometry in CEERS and PRIMER retaining only the available NIRCam wide filters (F090W, F115W, F150W, F200W, F277W, F356W, and F444W).F090W is not available in CEERS, which thus constitutes a more extreme test of the coverage of N U V at z > 3.In Figure B.4 we show the rest-frame N U V , U , V , and J flux densities in the two fitting runs, also removing the effect    2014).The color intensity scales as the density of points.Red filled and empty circles mark U V J-selected QGs from our sample with a counterpart at < 0. 2 and < 0. 5 in 3D-HST, respectively.From the left to right: photometric redshifts; stellar masses; total photometry in HST /F160W (in our analysis derived from the reference 0. 5 aperture); photometry in the same band in a common 0. 7 diameter aperture.The median offsets from the one-to-one relation (dotted lines) are shown, when applicable. of z phot .The results are fully consistent -and more so when F090W is included, as in the case of the test on PRIMER.

B.3. Sub-millimetric coverage and spectroscopically confirmed objects
We cross-checked our list of candidate quiescent objects with available catalogs of sub-millimetric surveys in CEERS (450 and 850 µm down to σ 450 = 1.2 and σ 850 = 0.2 mJy beam −1 with Scuba-2 in the deep tier of the S2CLS survey, Zavala et al. 2017; σ 850 = 0.2 mJy beam −1 over the full survey, Geach et al. 2017), PRIMER (870 µm with ALMA from the AS2UDS survey targeting Scuba-2 sub-millimeter galaxies from the S2CLS survey, Geach et al. 2017, and 2017) covers approximately 45% of our final samples in CEERS.The ALCS coverage of SMACS0723 is of ∼ 3 arcmin 2 centered on the cluster.AS2UDS and the ALMA archival observations are pointed and covered an area ∼ 600× smaller than the parent S2CLS survey in the UDS field (0.96 deg 2 , Stach et al. 2019).As mentioned in Section 3.4, we retrieve one ∼ 5σ-detection at 850 µm from Scuba-2 at a 0. 9 distance from a candidate U V J quiescent galaxy at z = 3.54 in CEERS (S2CLS-EGS-850.063 in Zavala et al. 2017, #9329 in our catalog).This candidate is selected by virtue of its uncertainty on the V − J color (0.4 mag) and the introduction of a padded box, while it is not picked by the N U V U V J criterion.However, several other possible optical/near-IR counterparts fall within the Scuba-2 beam (S.Gillman et al. in preparation), making the physical association inconclusive.Moreover, we matched our candidates with a compilation of spectroscopically confirmed galaxies from the literature.Despite the scarcity of these spectroscopic samples, we retrieve all sources in CEERS in both our selections from Schreiber et al. (2018b) and with fully consistent z phot ((z spec , z phot ): EGS-18996: (3.239, 3.12 +0.09 −0.05 ); EGS-40032: (3.219, 3.35 +0.09 −0.11 ; EGS-31322: (∼ 3.434, 3.54 +0.09 −0.10 )).We do not find any further matches with spectroscopically confirmed objects at any redshifts in our archive of Keck/MOSFIRE observations (G.Brammer et al. in preparation, Valentino et al. 2022) nor in the 3D-HST survey (Skelton et al. 2014;Momcheva et al. 2016).2022b) for a sample of 17 candidate QGs identified in CEERS by virtue of their low sSFR < 0.2/t obs , where t obs is the age of the Universe at the redshift of the galaxy.As mentioned in Section 3.4, there is an excellent overlap between our extended U V J selection and that in Carnall et al. (2022b), especially for their "robust" sample.Sources #9844, 4921 in our catalog (78374, 76507 in Carnall et al. 2022b) are excluded by virtue of their blue colors, while #9131 (92564) has a large uncertainty on V − J (σ V−J = 0.62 mag).The overlap is less extended when imposing P Q,50% ≥ 0.1.Sources below this threshold are either at the bluest (#9844, 4921) or reddest end of the color distribution (e.g., #7432, 8556 = 40015, 42128), the latter being mainly occupied by dusty SFGs.We remark the fact that our photometry is extracted in 0. 5 apertures and, thus, traces the properties of the central regions of galaxies.In presence of strong color gradients, as suggested by the RGB images of some of our candidates, photometry in larger apertures or based on surface brightness modeling across bands can drive to different results (e.g., #7432; see also Giménez-Arteaga Our photometry and SED modeling place this object at z phot = 3.58 +0.60 −0.24 and assigns it a M = 3.0 +1.4 −0.8 × 10 10 M and P Q,50% 0.1.
We highlight the fact that the comparison with both these works is partially affected by the different JWST zeropoint photometric calibration, an element in constant evolution to date.

C. STELLAR MASS LIMITS
Different mass limits could be a concern to draw comparison among fields with uneven photometric coverage and depth.In Figure C.6, we show that our compar- Open and red circles indicate the full and "robust" samples of candidates identified by the authors in the CEERS field.We did not apply any correction to homogenize our Chabrier (2003) IMF with that of Kroupa (2001).
ison is robust at the masses considered here.While a full-fledged analysis of the galaxy populations in our catalogs and stellar mass functions is deferred to future work, we show that the stellar mass cuts adopted in Section 4 are well above the threshold where completeness is expected to drop.For reference, we show our visually inspected U V J sample, more massive that the limit set by the 90% percentile of the M distribution in the 3 < z < 6.5 redshift interval in every field.This holds also for the robust N U V U V J-selected sample.
Figure C.7 shows the M distributions of the U V J and N U V − U, V − J selected sample.As noted in Section 3.3, candidates with P Q, 50% ≥ 0.7 maximally overlap with traditionally selected U V J galaxies.Both selection criteria tend to pick the most massive objects.Lower P Q, 50% thresholds also select bluer and less massive objects that might have recently quenched.We also note that our "strict" and "padded" subsamples of the overall robust, but looser visually inspected U V J pool of candidates do not introduce any immediately evident bias in mass.

D. LITERATURE COMPILATION
Table 4 includes complementary information on the comoving number densities of massive quiescent galaxies at 3 z 4 reported in Figure 5.The estimates are homogenized in terms of redshift interval and above the same mass limit of log(M /M ) 10.6 to the best of our knowledge.A lesser 0.1 dex difference in the mass integration limit remains for the estimates in Schreiber et al. (2018b) and Cecchi et al. (2019), while the redshift interval considered in Straatman et al. (2014) is 3.4 ≤ z < 4.2.The areas covered by the observations are those quoted in the original papers.There the reader can find further details about the exact selection techniques and their refinements, while we report the primary criterion in the table.The solid color lines mark the 90% percentile of the M distribution in redshift bins of equal number of points (i.e., 90% of the galaxies lie above these lines in each bin).For reference, the dotted lines indicate the percentile of the M distribution for U V J-selected red galaxies in the catalog at any redshifts (i.e., 90% of the U V J-selected QGs at any redshifts lie above these lines in each bin).We show the 90% percentile of the distribution of galaxies in the CEERS catalog in each panel (dashed blue line).A direct comparison of the 90% percentiles is shown in the bottom right panel.For reference, the red circles show the location of our visually-inspected U V J-selected sample of quiescent candidates at 3 < z < 6.5.

Figure 2 .
Figure 2. N U V − U , V − J rest-frame color diagrams for our combined sample of galaxies in JWST fields binned in redshift as labeled.Filled circles indicate our robust N U V U V J-selected quiescent sample (P Q,50% ≥ 0.1) and their 1σ uncertainties, color coded according to the nominal probability of being quiescent PQ for display purposes.The symbol size scales proportionally to the stellar mass as labeled.Thicker black circles show the sources with a robust or uncertain zspec in Schreiber et al. (2018b) falling in the portion of the CEERS field considered here.Gray points indicate the rest of the sample at those redshifts (Section 3).The color intensity scales as the density of points.The black arrow shows the effect of reddening for AV = 1.

Figure 4 .
Figure 4. Number densities of N U V − U , U − V , V − J-selected galaxies.The red and blue areas mark the p(z) of individual quiescent candidates at M ≥ 10 9.5 M with P Q,50% ≥ 0.7 and 0.1 ≤ P Q,50% < 0.7, respectively.The p(z) are normalized by their area ( z p(z) dz = 1).The sky coverage of each field and the comoving number densities per redshift bin in units of 10 −5 Mpc −3 are reported as in Figure3.The first and second rows indicate n in the 10 9.5 ≤ M < 10 10.6 M and M ≥ 10 10.6 M bins, respectively.
comoving number densities are expressed in units of 10 −5 Mpc −3 and computed over an area of 145.1 arcmin 2 .The uncertainties reflect the Poissonian 1σ confidence interval.

Figure 5 .
Figure 5. Comoving number densities of massive quiescent galaxies in the literature.The values have been homogenized in terms of redshift interval (3 z 4) and lower mass cut (log(M /M ) 10.6, similar IMF) to the largest possible extent.The uncertainties do not include the contribution of cosmic variance.The estimates are reported in Table4in Appendix D, along with complementary information.

Figure 6 .
Figure 6.Robust z > 4.5 quiescent candidates.Left column: spectral energy distributions.Black squares and blue filled circles indicate the observed and best-fit photometry of each source.Lighter gray squares mark observed flux densities with SNR < 3. Blue solid lines and shaded areas show the best-fit eazy-py models and their uncertainties.Central column: probability distribution functions of photometric redshifts z phot with eazy-py.The value of P Q, 50% is reported.Right column: SW and LW three color images of the candidates.The cutouts have sizes of 5 × 5 .

Figure A. 1 .
Figure A.1.Observed NIRCam F444W magnitudes as a function of the photometric redshifts.Gray points indicate sources in each field as label.The color intensity scales as the density of points.Red circles show our U V J-selected sample of quiescent candidates at 3 < z < 6.5 after the visual inspection.The color lines mark the 5σ depths in 0. 5 diameter apertures in F444W.For reference, we show the depth for the CEERS field in each panel (dashed blue line).A direct comparison of the depth is shown in the bottom right panel.MIRI.Here we considered the area covered with NIR-Cam in the UDS field, critically overlapping with the HST deep imaging from the Cosmic Assembly Nearinfrared Deep Extragalactic Legacy Survey (CANDELS, Grogin et al. 2011).A.4. North Ecliptic Pole The North Ecliptic Pole (NEP) Time-Domain Field (TDF) is being observed as part of the Guaranteed-Time Observations (GTO, program 2738, PI: R. Windhorst).

Figure
Figure A.2. JWST coverage maps.For every field, we show the footprint of each JWST filter colored as labeled.The red shaded area indicates the overlap of our selection filters (F150W, F200W, F277W, F356W, and F444W).

Figure B. 3 .
Figure B.3.Comparison with 3D-HST.Gray points indicate sources in common (maximal separation < 0. 5) between our catalogs in CEERS (top row) and PRIMER (bottom row) and those in 3D-HST from Skelton et al. (2014).The color intensity scales as the density of points.Red filled and empty circles mark U V J-selected QGs from our sample with a counterpart at < 0. 2 and < 0. 5 in 3D-HST, respectively.From the left to right: photometric redshifts; stellar masses; total photometry in HST /F160W (in our analysis derived from the reference 0. 5 aperture); photometry in the same band in a common 0. 7 diameter aperture.The median offsets from the one-to-one relation (dotted lines) are shown, when applicable.

Figure B. 4 .
Figure B.4.Effect of HST filter coverage.From left to right: rest-frame N U V , U , V , and J flux densities from SED modeling with eazy-py in the CEERS field.Quantities on the X-axes are computed with the full filter set, including HST bands; those on the Y-axes are from the run with NIRCam wide filters only.The flux densities are computed at their respective best-fit z phot in the top row and at fixed z phot (full filter set) in the bottom one to remove the effect of different redshifts.Gray dots indicate the whole galaxy sample.Blue dots mark bright objects with F200W < 28 mag in the redshift interval 3 ≤ z ≤ 6.5 of interest, as labeled.Large filled circles indicate our selected QGs at z > 3 using the N U V U V J diagram with their uncertainties, color coded according to the nominal PQ in the run with the full filter set.B.4.Comparison with JWST-selected photometric quiescent candidates in the literature Figure B.5 shows the comparison between our F200W magnitudes and SED modeling results with eazy-py and those from Carnall et al. (2022b) for a sample of 17 candidate QGs identified in CEERS by virtue of their low sSFR < 0.2/t obs , where t obs is the age of the Universe at the redshift of the galaxy.As mentioned in Section 3.4, there is an excellent overlap between our extended U V J selection and that inCarnall et al. (2022b), especially for their "robust" sample.Sources #9844, 4921 in our catalog(78374, 76507 in Carnall  et al. 2022b) are excluded by virtue of their blue colors, while #9131 (92564) has a large uncertainty on V − J (σ V−J = 0.62 mag).The overlap is less extended when imposing P Q,50% ≥ 0.1.Sources below this threshold are either at the bluest(#9844, 4921)  or reddest end of the color distribution(e.g., #7432, 8556 = 40015, 42128), the latter being mainly occupied by dusty SFGs.We remark the fact that our photometry is extracted in 0. 5 apertures and, thus, traces the properties of the central regions of galaxies.In presence of strong color gradients, as suggested by the RGB images of some of our candidates, photometry in larger apertures or based on surface brightness modeling across bands can drive to different results (e.g., #7432; see also Giménez-Arteaga et al. 2022).Despite this, we find an overall agreement in z phot and M (Figure B.5).If any, our z phot seem to be systematically lower and M larger than those derived by Carnall et al. (2022b) (Kocevski et al. 2022 also report lower redshift estimates).However, these offsets are in the realm of typical statistical and systematic uncertainties that different codes ran with a variety of parameters can produce.In addition, our selections do not retrieve the dusty candidate QG at z phot ∼ 5.4 in SMACS0723 presented in Rodighiero et al. (2023, ID#2=KLAMA, #1536 in our catalog; R.A. = 110.70257564,Dec. = −73.48472291 in the Gaia DR3 astrometric reference).

Figure B. 5 .
Figure B.5.Comparison with JWST -selected QG at z > 3 in CEERS.Photometric redshift z phot , M , and NIRCam F200W in our analysis (Y-axis) and in that presented in Carnall et al. (2022b) (X-axes).Open and red circles indicate the full and "robust" samples of candidates identified by the authors in the CEERS field.We did not apply any correction to homogenize ourChabrier (2003) IMF with that ofKroupa (2001).

Figure C. 6 .
Figure C.6.Stellar masses as a function of the photometric redshifts.Gray points indicate sources in each field as label.The color intensity scales as the density of points.The solid color lines mark the 90% percentile of the M distribution in redshift bins of equal number of points (i.e., 90% of the galaxies lie above these lines in each bin).For reference, the dotted lines indicate the percentile of the M distribution for U V J-selected red galaxies in the catalog at any redshifts (i.e., 90% of the U V J-selected QGs at any redshifts lie above these lines in each bin).We show the 90% percentile of the distribution of galaxies in the CEERS catalog in each panel (dashed blue line).A direct comparison of the 90% percentiles is shown in the bottom right panel.For reference, the red circles show the location of our visually-inspected U V J-selected sample of quiescent candidates at 3 < z < 6.5.

Figure C. 7 .
Figure C.7. Stellar mass distributions of the selected samples.The solid lines mark the probability density functions of M for the U V J and N U V U V J final samples colored as labeled.The curves are smoothed with a Gaussian Kernel Density Estimator and normalized to an area of 1. Top row: Redshift interval 3 < z < 4. Bottom row: 4 < z < 5.Only # 185 in PRIMER is selected with P Q, 50% ≥ 0.7 at z > 4. The red arrow marks its stellar mass estimate in the bottom right panel.

Table 1 .
Properties of the observed fields with JWST /NIRCam observations.

Table 2 .
Comoving number densities of quiescent galaxies in this work.

Table 3 .
Filter coverage in each field.
. Figure A.1 shows the depths in F444W within the apertures used for the photometric extraction.Similar plots for F150W, F200W, F277W, and F356W are available in Figure Set A1.The depths are reported in Table 1.The minimum overlap of the NIRCam bands imposed for the selection and corresponding to the areas in Table 3 is shown in Figure A.2. HST legacy field with several JWST instruments for imaging and, in the future, spectroscopy

Table 4 .
Comoving number densities of quiescent galaxies at 3 < z < 4 in the literature.