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

Francesco Valentino, Gabriel Brammer, Katriona M. L. Gould, Vasily Kokorev, Seiji Fujimoto, Christian Kragh Jespersen, Aswin P. Vijayan, John R. Weaver, Kei Ito, Masayuki Tanaka, et al.

► To cite this version:

Francesco Valentino, Gabriel Brammer, Katriona M. L. Gould, Vasily Kokorev, Seiji Fujimoto, et al. An Atlas of Color-selected Quiescent Galaxies at $z > 3$ in Public JWST Fields. The Astrophysical Journal, 2023, 947, 10.3847/1538-4357/acbefa. insu-04479092

HAL Id: insu-04479092

https://insu.hal.science/insu-04479092

Submitted on 4 Mar 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

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

CrossMark

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

Francesco Valentino<sup>1,2,3</sup>, Gabriel Brammer<sup>1,2</sup>, Katriona M. L. Gould<sup>1,2</sup>, Vasily Kokorev<sup>1,4</sup>, Seiji Fujimoto<sup>1,5</sup>, Christian Kragh Jespersen<sup>6</sup>, Aswin P. Vijayan<sup>1,7</sup>, John R. Weaver<sup>8</sup>, Kei Ito<sup>9</sup>, Masayuki Tanaka<sup>10,11</sup>, Olivier Ilbert<sup>12</sup>, Georgios E. Magdis<sup>1,2,7</sup>, Katherine E. Whitaker<sup>1,8</sup>, Andreas L. Faisst<sup>13</sup>, Anna Gallazzi<sup>14</sup>, Steven Gillman<sup>1,7</sup>, Clara Giménez-Arteaga<sup>1,2</sup>, Carlos Gómez-Guijarro<sup>15</sup>, Mariko Kubo<sup>16,17</sup>, Kasper E. Heintz<sup>1,2</sup>, Michaela Hirschmann<sup>18,19</sup>, Pascal Oesch<sup>1,2,20</sup>, Masato Onodera<sup>11,21</sup>, Francesca Rizzo<sup>1,2</sup>, Minju Lee<sup>1,7</sup>, Victoria Strait<sup>1,2</sup>, and Sune Toft<sup>1,2</sup>

<sup>1</sup> Cosmic Dawn Center (DAWN), Denmark; francesco.valentino@eso.org

<sup>2</sup> Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark

<sup>3</sup> European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei Munchen, Germany

<sup>4</sup> Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700AV Groningen, The Netherlands

<sup>5</sup> Department of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA

<sup>6</sup> Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

<sup>7</sup> DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800, Kgs. Lyngby, Denmark

<sup>8</sup> Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA

<sup>9</sup> Department of Astronomy, School of Science, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan

<sup>10</sup> National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo, 181-8588, Japan

<sup>11</sup> Department of Astronomical Science, The Graduate University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo, 181-8588, Japan

<sup>12</sup> Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France

<sup>13</sup> Caltech/IPAC, MS314-6, 1200 E. California Blvd. Pasadena, CA 91125, USA

<sup>14</sup> INAF—Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5, I-50125 Firenze, Italy

<sup>15</sup> Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191, Gif-sur-Yvette, France

<sup>16</sup> Astronomical Institute, Tohoku University, Aramaki, Aoba-ku, Sendai, Miyagi 980-8578, Japan

<sup>17</sup> Research Center for Space and Cosmic Evolution, Ehime University, 2-5 Bunkyo-cho, Matsuyama, Ehime 790-8577, Japan

<sup>18</sup> Institute for Physics, Laboratory for Galaxy Evolution and Spectral Modelling, Ecole Polytechnique Federale de Lausanne, Observatoire de Sauverny, Chemin Pegasi 51, 1290 Versoix, Switzerland

<sup>19</sup> INAF—Osservatorio Astrometrico di Trieste, Via Tiepolo 11, I-34131 Trieste, Italy

<sup>20</sup> Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland

<sup>21</sup> Subaru Telescope, National Astronomical Observatory of Japan, National Institutes of Natural Sciences (NINS), 650 North A'ohoku Place, Hilo, HI 96720, USA
Received 2023 January 3; revised 2023 February 9; accepted 2023 February 9; published 2023 April 14

## Abstract

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 arcmin². 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_{\star} \ge 10^{10.6} M_{\odot}$) 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.

*Unified Astronomy Thesaurus concepts: Galaxy evolution (594); High-redshift galaxies (734); Galaxy quenching (2040); Quenched galaxies (2016); Post-starburst galaxies (2176); Surveys (1671)*

*Supporting material: figure set, tar.gz file*

## 1. Introduction

Over the past few years, the existence of a population of quenched and quiescent galaxies (QGs) at redshifts $z \sim 3-4$ (e.g., Fontana et al. 2009; Spitler et al. 2014; Straatman et al. 2014) has been finally corroborated by the long sought after spectroscopic confirmations (Glazebrook et al. 2017;

Schreiber et al. 2018a, 2018b; Tanaka et al. 2019; D'Eugenio et al. 2020, 2021; Forrest et al. 2020a, 2020b; Valentino et al. 2020; 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 (SFRs), also supported with long-wavelength observations (Santini et al. 2019, 2021; Suzuki et al. 2022); emission from active galactic nuclei potentially pointing at a coevolution with or feedback from their central supermassive black holes (Marsan et al. 2015, 2017; Ito et al. 2022; Kubo et al. 2022); stellar velocity

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

Original content from this work may be used under the terms
of the Creative Commons Attribution 4.0 licence. Any further
distribution of this work must maintain attribution to the author(s) and the title
of the work, journal citation and DOI.
dispersions ($\sigma_*$) and dynamical masses (Tanaka et al. 2019; Saracco et al. 2020) with possible implications on their initial mass function (IMF; 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 submillimeter galaxies with depletion timescales of $\tau \lesssim 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 \sim 2$ QGs (Cimatti et al. 2008; Toft et al. 2014; Akhshik et al. 2023). 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” submillimeter-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 (Ciesla et al. 2016; Carnall et al. 2018, 2019; Schreiber et al. 2018a; Leja et al. 2019a; Iyer et al. 2019; K. Gould et al. 2023, 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; Long et al. 2022; Manning 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 \gtrsim 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 $\ge 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 megayears) as in “post-starburst” galaxies rather than being prototypical “red and dead” objects (Schreiber et al. 2018a; D’Eugenio et al. 2020; Forrest et al. 2020b; Lustig et al. 2021; Gould et al. 2023; Marsan et al. 2022), even if examples of older populations are available (Glazebrook et al. 2017; M. Tanaka et al. 2023, 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 simultaneous 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 \sim 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 $\sigma_*$ 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 have prevented us from unambiguously confirming whether QGs exist at $z > 4$ (see Merlin et al. 2019; Carnall et al. 2020; Mawatari et al. 2020), 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 (see Carnall et al. 2023c). 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 multi-wavelength 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 $\sim 145$ arcmin$^2$ at $z \sim 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. 2018a; Merlin et al. 2019). In addition, 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 $\Lambda$CDM cosmology with $\Omega_m = 0.3$, $\Omega_\Lambda = 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.$^{22}$

## 2. 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. 2023, in preparation). The approach is similar to that in Labbe et al. (2022) and Bradley et al. (2022), here including the recently updated zero-points.

<footnote>$^{22}$ Supplemental data are made available on Zenodo doi:10.5281/zenodo.7614908, and the JWST and HST image mosaics created with grizli are stored on the Electronic Research Data Archive at the University of Copenhagen (Brammer 2023).</footnote>
<table><caption>**Table 1**

Properties of the Observed Fields with JWST/NIRCam Observations</caption><thead><tr><th>Field</th><th>R.A.<br>(deg)</th><th>Decl.<br>(deg)</th><th>Area<br>(arcmin<sup>2</sup>)</th><th>NIRCam Depths<br>(mag)</th><th>HST</th></tr></thead><tbody><tr><td>CEERS</td><td>214.88598</td><td>52.89500</td><td>34.7</td><td>28.5 / 28.8 / 28.8 / 28.8 / 28.3</td><td>Yes</td></tr><tr><td>Stephan's Quintet</td><td>339.00057</td><td>33.95996</td><td>35.0<sup>a</sup></td><td>27.5 / 27.6 / 28.0 / 28.1 / 27.7</td><td>No</td></tr><tr><td>PRIMER</td><td>34.37792</td><td>-5.14717</td><td>21.9</td><td>27.5 / 27.7 / 27.9 / 27.9 / 27.4</td><td>Yes</td></tr><tr><td>NEP</td><td>260.73773</td><td>65.78167</td><td>9.7</td><td>28.5 / 28.6 / 28.9 / 28.9 / 28.3</td><td>Yes</td></tr><tr><td>J1235</td><td>188.96741</td><td>4.92465</td><td>9.0</td><td>28.4 / 29.1 / 29.3 / 29.3 / 28.4</td><td>No</td></tr><tr><td>GLASS</td><td>3.50145</td><td>-30.33612</td><td>8.5</td><td>28.8 / 29.0 / 29.1 / 29.1 / 29.4</td><td>Yes</td></tr><tr><td>Sunrise</td><td>24.34743</td><td>-8.43215</td><td>7.3<sup>b</sup></td><td>28.1 / 28.3 / 28.4 / 28.4 / 28.0</td><td>Yes</td></tr><tr><td>SMACS 0723</td><td>110.75478</td><td>-73.46788</td><td>6.5<sup>b</sup></td><td>28.8 / 29.0 / 29.2 / 29.2 / 28.8</td><td>Yes</td></tr><tr><td>SGAS 1723</td><td>260.91450</td><td>34.19371</td><td>5.3</td><td>25.8 / 25.9 / 26.6 / 26.7 / 26.6</td><td>Yes</td></tr><tr><td>SPT 0418</td><td>64.66113</td><td>-47.87526</td><td>5.0</td><td>26.6 / 27.1 / 27.8 / 27.3 / 27.1</td><td>No</td></tr><tr><td>SPT 2147<sup>c</sup></td><td>326.82917</td><td>-50.59632</td><td>2.3</td><td>— / 27.7 / 27.4 / 27.7 / 26.9</td><td>Yes</td></tr></tbody></table>

Notes. 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 \sim 3-5$ (Section 4).

c No F150W coverage.

## 2.1. Reduction

We homogeneously process the publicly available JWST imaging obtained with the NIRCam, NIRISS, and MIRI instruments in 11 fields targeted during the first 3 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 zero-points relative to *jwst_0942.pmap*, including detector variations (Brammer 2022). 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.pmap* calibration 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 NIRCam 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), coadd, and finally drizzle them (Fruchter & Hook 2002) to a 0.′′02 pixel scale for the Short Wavelength (SW) NIRCam bands and to 0.′′04 for all the remaining JWST and Hubble Space Telescope (HST) filters. We provide further details about the individual fields in Appendix A.

## 2.2. Extraction

We extract sources using a detection image produced by combining all of the “wide” (W) NIRCam 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).<fnref n="23" /> The aperture correction is computed on the LW

detection image and applied to all bands. The depths in the reference 0.′′5 apertures in the five NIRCam bands that we require to select candidate QGs (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 Figure A1 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 (PSFs).<fnref n="24" />

## 2.3. Modeling of the Spectral Energy Distribution

We utilize EAZY-PY<fnref n="25" /> (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 zero-point 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 maximum flexibility. This set of templates covers a large interval in ages, dust attenuation, and lognormal SFHs—spanning the whole UVJ rest-frame color diagram. More specifically, the *corr_sfhz_13* subset of models within EAZY contains redshift-dependent SFHs, which, at a given redshift, exclude 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. (2023a)—is also included to allow for an extra degree of freedom in photometric solutions of distant objects, but it is not accounted for in the stellar mass calculation.<fnref n="26" /> The templates are created adopting a Chabrier (2003) IMF and applying a Kriek & Conroy (2013) dust attenuation law (dust index

<footnote n="23">²³ Photometric measurements in different apertures are available in the online catalogs (see footnote 22).</footnote>

<footnote n="24">²⁴ https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/psf</footnote>

<footnote n="25">²⁵ https://github.com/gbrammer/eazy-py</footnote>

<footnote n="26">²⁶ The parameters attached to this template are uncertain (Carnall et al. 2023a; Giménez-Arteaga et al. 2022). Solutions dominated by this template should be treated with caution. This is not the case for galaxies in our final samples.</footnote>
$\delta = -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_\odot) \in [-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. Thus, we 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.<fnref n="27" /> In addition, 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 SMACS 0723 ($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 of photometric redshifts, we obtain a good agreement with spectroscopic determinations from archival observations when available (G. Brammer et al. 2023, in preparation). For reference, in the CEERS field we estimate a $\sigma_{NMAD} = 0.0268$ (0.0187) for the spectroscopic sample (clipping catastrophic outliers).<fnref n="28" /> Stellar masses are consistent with independent estimates obtained with finer grid-based codes (see Figure B1 for a comparison with 3D-HST). SFRs 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.

## 2.4. 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 16th–84th percentiles of 100 fits drawn from the best-fit template error function. We compute the rest-frame magnitudes in the GALEX near-UV (NUV) band ($\lambda = 2800 \AA$), the *U* and *V* Johnson filters defined as in Maíz Apellániz (2006), and the *J* Two Micron All Sky Survey 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; Appendix A in Gould et al. 2023). This allows for a color determination that relieves the dependency on the adopted models, while using the whole photometric information.

## 3. Sample Selection

Before selecting QG 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 $\chi^2/N_{filt} \le 8$, where $N_{filt} \ge 6$ is the number of available filters. The latter includes NIRCam wide bands at 1.5, 2.0, 2.7, 3.5, and 4.4 $\mu$m in every field with the exception of SPT 2147, 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 $\mu$m allows for robust determinations of stellar masses. Then, we apply a loose cut at $>5$ on the quadratic sum of the signal-to-noise ratio 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\%}$) / ($2p(z)_{50\%}$) $< 0.3$, where $p(z)_{1\%}$ indicates the $i$th percentile of $p(z)$). Finally, we apply a cut in redshift at $3 \le z \le 6.5$ with a buffer of $dz = 0.1$ and accounting for the uncertainty on the best-fit solution ($p(z)_{84\%} \ge 3 - dz$ and $p(z)_{16\%} \le 6.5 + dz$). To pick quiescent objects, we opt for a rest-frame color selection following a dual approach.

### 3.1. UVJ Color Diagram

On the one hand, we select objects in the classical UVJ diagram. This allows us to directly compare our results with a large body of literature that has accumulated over the past 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\sigma$ uncertainties on the colors consistent with the selection box as long as $\sigma_{color} < 0.5$ mag. We then visually inspect the images and the spectral energy distribution (SED) fits of 251 candidate QGs. 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 owing 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 an 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\sigma$ 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 UVJ-selected samples.

### 3.2. NUV - U, V - J Color Diagram

In parallel, we follow the novel method described in Gould et al. (2023; see also Antwi-Danso et al. 2023 for an alternative approach introducing a synthetic band). The authors incorporate the NUV magnitude in their selection and model the galaxy distribution in the $(NUV - 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 NUV 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

<footnote n="27">²⁷ https://github.com/gbrammer/eazy-photoz/tree/master/templates/sfhz</footnote>

<footnote n="28">²⁸ The redshift distributions and the comparison with spectroscopic estimates (e.g., $\Delta z/(1+z)$ as a function of $z$) for each field are bundled with the mosaics (see footnote 22).</footnote>
<figure><img src="image_4.png" /><figcaption>**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 *UVJ* quiescent sample and their $1\sigma$ 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 $A_V = 1$.</figcaption></figure>

epoch of quenching of the first galaxies (D'Eugenio et al. 2020; Forrest et al. 2020b). Moreover, the GMM allows us to fully account for the blurred separation between star-forming galaxies and QGs 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. 2022b) assuming $5 \times$ 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 16th–84th percentiles of the distribution as our reference $P_{Q,50\%}^{29}$ and its $1\sigma$ 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\%} \ge 0.1$ to select candidate QGs, 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 *UVJ*-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 \le P_{Q,50\%} < 0.7$ and 18/20 (~90%) truly passive galaxy candidates with $P_{Q,50\%} \ge 0.7$. Their location in the projected NUV–*U*, *V*–*J* plane is shown in Figure 2.

### 3.3. Overlap between Selections

As noted in Gould et al. (2023), a selection in the NUV *UVJ* arguably outperforms the classical *UVJ* in selecting quiescent (passive and recently quenched or post-starburst) 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 *UVJ* 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 NUV *UVJ* and *UVJ*-selected objects, respectively, comparable to 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\%} \ge 0.7$ and $0.1 \le P_{Q,50\%} < 0.7$ are part of the *UVJ* sample. Moreover, 16/18 objects with $P_{Q,50\%} \ge 0.7$ fall within the standard *UVJ* 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 that resemble more classically QGs owing to how the model is trained (Figure 2). Therefore, it is expected that the *UVJ*-selected sample has a smaller overlap with galaxies that have lower $P_{Q,50\%}$ values and that the galaxies that it does not select are those that are recently quenched. The overlap is reflected also in the $M_*$ distributions of the selected samples (see Figure C2). Lower $P_{Q,50\%}$ values are associated with bluer, more recently quenched, but also lower-mass objects, otherwise missed by *UVJ* selections.

### 3.4. Sanity Checks on the Sample

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

#### 3.4.1. 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 B1). 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_\odot)$ (<0.2 dex) make our results fully consistent between different SED modeling codes (Appendix B.1).

#### 3.4.2. Availability of HST Imaging

We also test the impact on our sample selection of HST filters, which increase the sampling of the rest-frame 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 $\mu$m, mimicking the situation for Stephan's Quintet, where no HST imaging is at disposal. We obtain fully consistent $z_{phot}$, $M_*$, and

<footnote n="29">29 For clarity, we stress that, in this notation, “50%” refers to the percentile of the bootstrapped $P_Q$ and not to a probability of 50% to be quiescent.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** NUV–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 NUVUVJ-selected quiescent sample ($P_{Q,50\%} \ge 0.1$) and their $1\sigma$ uncertainties, color-coded according to the nominal probability of being quiescent $P_Q$ 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 $z_{\text{spec}}$ in Schreiber et al. (2018a) 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 $A_V = 1$.</figcaption></figure>

rest-frame color estimates within the uncertainties, especially when removing the effect of $z_{\text{phot}}$ from the calculation and focusing on the $3 \le z \le 6.5$ interval of interest (Figure B2). This holds also when F090W, probing wavelengths shorter than NUV at the lower end of the redshift interval that we explore, is not available, as in the case of CEERS. We also reapply the NUVUVJ selection, including bootstrapping, and obtain consistent samples taking into account the uncertainties.

### 3.4.3. 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 submillimetric observations in CEERS (450 and 850 µm with SCUBA-2; Geach et al. 2017; Zavala et al. 2017), PRIMER (870 µm with the Atacama Large Millimeter/ submillimeter Array (ALMA) from the AS2UDS survey; Dudzevičiūtė et al. 2020, and Cheng et al. 2023), and SMACS 0723 (1.1 mm with ALMA from the ALCS Survey; Kokorev et al. 2022; Fujimoto et al. 2023). 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 submillimeter galaxy candidate possibly associated with an overdensity (S. Gillman et al. 2023, in preparation). Removing it from our sample does not change the results of this work. The absence of long-wavelength 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).

### 3.4.4. Spectroscopically Confirmed and Alternative JWST-based Photometric Samples

We correctly identify and select the spectroscopically confirmed QGs in Schreiber et al. (2018a) at consistent $z_{\text{phot}}$. Our relaxed UVJ cuts and $P_{Q,50\%} \ge 0.1$ recover 14/17 and 11/17 candidate QGs selected via SED modeling and a specific SFR (sSFR) threshold in Carnall et al. (2023b), respectively. When considering only their “robust” sample, there is a 89% and 78% overlap with our UVJ and $P_{Q,50\%}$ selections, respectively. We calculate lower but generally consistent $z_{\text{phot}}$ (Figure B3; see also Kocevski et al. 2023). Minor systematic offsets in the $M_*$ and F200W magnitudes are present (Figure B3).

## 4. Number Densities

For each field, we compute the comoving number density $n$ of candidate QGs in three redshift bins ($z \in [3, 4], [4, 5], \text{ 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_{\text{phot}}$ within their $p(z)$ and compute the median and 16th–84th 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 SMACS 0723 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, 2021) using the available HST and JWST data (Harikane et al. 2023; Welch et al. 2022b). The effect of lensing on the effective area varies negligibly in the redshift interval and luminosity regime spanned by our samples of candidate QGs. The linear magnification for the only QG candidate in proximity of SMACS 0723 is ~1.7. For two candidates in WHL 0137 (Sunrise) with $\log(M_*/M_\odot) < 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 UVJ- and NUVUVJ-selected QGs in each of the 11 fields that we consider. Combined estimates based on the aggregate area of 145.1 arcmin² are also presented and reported in Table 2. In each panel, we report the number densities in stellar mass bins of $10^{9.5} M_\odot \le M_* < 10^{10.6} M_\odot$ and $M_* \ge 10^{10.6} M_\odot$. The high mass threshold is chosen to directly compare these results with those in the literature (Table 4). Such a threshold also allows us to safely compare different fields. This is clear from Figure C1, showing the stellar mass limit in each field for the overall sample of galaxies and QGs.

For the UVJ-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
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Number densities of UVJ-selected galaxies in public JWST fields. The purple and green colors mark the $p(z)$ of individual robust “strict” UVJ quiescent candidates with $M_* \ge 10^{10.6} M_\odot$ and $10^{9.5} M_\odot \le M_* < 10^{10.6} M_\odot$, respectively. The $p(z)$ are normalized by their area ($\int_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, SMACS 0723), for which the masses should be intended as observed (not delensed). The comoving number densities in units of $10^{-5} \text{ 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_\odot \le M_* < 10^{10.6} M_\odot$ and $M_* \ge 10^{10.6} M_\odot$ bins, respectively.</figcaption></figure>

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 UVJ color box returns $2 \times$ and $1.2 \times$ 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 UVJ selection; Schreiber et al. 2018a; thus likely higher for the padded and the robust samples). The NUVUVJ-selected sample

$(P_{Q,50\%} \ge 0.1)$ is as numerous as the UVJ-selected one, despite the partial overlap between the two criteria (Section 3.3), similar to what is found in Gould et al. (2023).

## 4.1. 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_\odot)$ 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 (Chuang & Lin 2023; Jespersen et al. 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
<figure><img src="image_7.png" /><figcaption>Figure 4. Number densities of NUV - 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 Figure 3. The first and second rows indicate n in the 10^9.5 M⊙ ≤ M* < 10^10.6 M⊙ and M* ≥ 10^10.6 M⊙ bins, respectively.</figcaption></figure>

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. (2022a).
The total cosmic variance is then computed as

$$
\sigma_{\text{CV, total}} = \sqrt{\frac{1}{\sum_{\text{fields}} \sigma_{\text{CV, field}}^2}}, \quad (1)
$$

which assumes that all fields are independent.<fnref n="30" /> 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; Cecchi et al. 2019; Girelli
et al. 2019; Merlin et al. 2019; Shahidi et al. 2020; Casey et al.
2022; Gould et al. 2023; Long et al. 2022; Weaver et al. 2022a).
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
(Kroupa 2001; Chabrier 2003). We also add number density
estimates from EAGLE (Crain et al. 2015; Schaye et al. 2015),
Illustris (Vogelsberger et al. 2014), IllustrisTNG100, and
IllustrisTNG300 simulations (Nelson et al. 2019). We count

<footnote n="30">³⁰ Note that if all fields were the same shape and area, this formula reduces to the well-known $\sigma_{\text{CV, total}} = \sigma_{\text{CV, field}} / \sqrt{N}$.</footnote>
<table><caption>**Table 2**
Comoving Number Densities of Quiescent Galaxies in This Work</caption><thead><tr><th>Redshift</th><th>log(M<sub>*</sub>)<br>(M<sub>☉</sub>)</th><th>UVJ<br>Strict</th><th>UVJ<br>Padded</th><th>NUVUVJ<br>P<sub>Q,50%</sub></th><th>σ<sub>CV</sub><br>(%)</th></tr></thead><tbody><tr><td rowspan="2">3 &lt; z &lt; 4</td><td>[9.5, 10.6)</td><td>3.9<sup>+1.2</sup><sub>-0.9</sub></td><td>4.1<sup>+1.2</sup><sub>-0.9</sub></td><td>2.8<sup>+1.0</sup><sub>-0.8</sub></td><td>0.10</td></tr><tr><td>&gt;10.6</td><td>2.4<sup>+1.0</sup><sub>-0.7</sub></td><td>2.7<sup>+1.0</sup><sub>-0.8</sub></td><td>2.3<sup>+1.0</sup><sub>-0.7</sub></td><td>0.18</td></tr><tr><td rowspan="2">4 &lt; z &lt; 5</td><td>[9.5, 10.6)</td><td>0.6<sup>+0.7</sup><sub>-0.3</sub></td><td>1.0<sup>+0.8</sup><sub>-0.5</sub></td><td>1.0<sup>+0.8</sup><sub>-0.5</sub></td><td>0.16</td></tr><tr><td>&gt;10.6</td><td>0.7<sup>+0.7</sup><sub>-0.4</sub></td><td>0.9<sup>+0.8</sup><sub>-0.5</sub></td><td>0.9<sup>+0.7</sup><sub>-0.4</sub></td><td>0.30</td></tr><tr><td rowspan="2">5 &lt; z &lt; 6.5</td><td>[9.5, 10.6)</td><td>0.0<sup>+0.3</sup><sub>-0.0</sub></td><td>0.2<sup>+0.4</sup><sub>-0.2</sub></td><td>0.1<sup>+0.4</sup><sub>-0.1</sub></td><td>0.22</td></tr><tr><td>&gt;10.6</td><td>0.0<sup>+0.3</sup><sub>-0.0</sub></td><td>0.0<sup>+0.3</sup><sub>-0.0</sub></td><td>0.0<sup>+0.3</sup><sub>-0.0</sub></td><td>0.41</td></tr></tbody></table>

**Note.** The 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\sigma$ confidence interval. Upper limits are at $1\sigma$ using the same approach (Gehrels 1986). Statistical uncertainties are accounted for by integrating the $p(z)$ within the redshift intervals. The uncertainties due to cosmic variance are expressed as fractional $\sigma_{CV}$ deviations (Section 4.1). The selections are described in Section 3. The adopted threshold for the NUVUVJ selection is $P_{Q,50\%} \ge 0.1$.

simulated galaxies with sSFR $\le 10^{-10}$ yr$^{-1}$ within 4 × and 2 × the half-mass radius for EAGLE 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 $\sim 2.5 \times 10^{-5}$ Mpc$^{-3}$, consistent with some of the determinations with similar color or sSFR cuts (Schreiber et al. 2018a; Merlin et al. 2019). Our estimates are $\sim 2 \times$ larger than the most recent measurements in the largest contiguous survey among those considered, COSMOS (Weaver et al. 2022a), 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; Cecchi et al. 2019; Girelli et al. 2019). This is due to a combination of deeper and homogeneous measurements in the optical and near-IR over a twice as large effective area, now available in COSMOS2020 (Weaver et al. 2022b), more conservative and pure samples of QGs, the specific templates used in each work, and the integration of a 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 allow 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 QGs 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. 2023b). Nevertheless, in our most massive bin, the variations among different works are still dominated by systematics in the selection, modeling, and cosmic variance.

### 4.3. 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 $\sim 145$ arcmin$^2$, we find per-field value oscillations of a factor of $2-3 \times$ (Table 5). We ascribe these differences to cosmic

variance and to the fact that massive quiescent systems might already be signposts of distant overdensities and protoclusters—massive halos able to fast-track galaxy evolution. In Table 5, we report the fractional $1\sigma$ uncertainties due to cosmic variance in each field. Taken individually, an uncertainty of $\sim 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 overdense despite a similar sky coverage. CEERS displays the largest number densities of QGs $> 10^{10.6} M_\odot$ in our compilation, consistent with the estimates in Carnall et al. (2023b). 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. 2023b, respectively). The pair, possibly interacting, is surrounded by two more red sources with similar $z_{phot} \sim 3.18–3.54$ that fall in the visually vetted UVJ 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 \gtrsim 3$ (Kalita et al. 2021; Kubo et al. 2021; McConachie et al. 2022) or in close pairs with other massive objects (“Jekyll & Hyde” at $z = 3.717$; Glazebrook et al. 2017; Schreiber et al. 2018b; 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 UVJ-selected sample.

### 4.4. A Look at Lower Masses

Figures 3 and 4 show the number densities at lower stellar masses ($10^{9.5} M_\odot \le M_* < 10^{10.6} M_\odot$), 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 $\sim 10^8 M_\odot$ at $3 < z < 6.5$ (excluding SPT 2147; see Figure C1). The lower limit of $M_* = 10^{9.5} M_\odot$ 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 ($\sim 1.2–1.6 \times$) 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 QGs, roughly peaking and flattening or turning over at $\sim 10^{10.6} M_\odot$ at these redshifts (Weaver et al. 2022a) and revealing a steady buildup of lower-mass QGs (Santini et al. 2022). Promising low-mass QGs have already been confirmed with JWST/NIRISS at $z \sim 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.

### 4.5. 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 $\sim 2.5–4.5 \times$ smaller above and below $M_* = 10^{10.6} M_\odot$ than those in CEERS by Carnall et al. (2023b), who also perform a selection starting from JWST images. The difference shrinks to a factor of $\sim 2.5–1 \times$ when we consider only the same field and the “robust” sample in
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comoving number densities of massive QGs in the literature. The values have been homogenized in terms of redshift interval ($3 \lesssim z \lesssim 4$) and lower mass cut ($\log(M_*/M_\odot) \gtrsim 10.6$, similar IMF) to the largest possible extent. The uncertainties do not include the contribution of cosmic variance. The estimates are reported in Table 4 in Appendix D, along with complementary information.</figcaption></figure>

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_* \ge 10^{10.6} M_\odot$) are consistent with those in the COSMOS field by Weaver et al. (2022a), but $1.8 \times$ larger than what was 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_\odot$, but not homogenized among different works at this stage, given the impact of different depths at $z>4$), we retrieve similar $n$ to that in COSMOS ($(1.0 \pm 0.3) \times 10^{-5} Mpc^{-3}$ for $M_* > 10^{9.9} M_\odot$; Weaver et al. 2022a) and $1.6-2 \times$ larger than in large-field HST surveys such as CANDELS ($\sim 7.9 \times 10^{-6} Mpc^{-3}$ for the “complete” sample at $M_* > 5 \times 10^9 M_\odot$; 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 five most robust sources falling in either the “strict” or “padded” UVJ selections are shown in Figure 6. We do not find reliable candidates at $z \gtrsim 5.2$, which signposts the earliest epoch of appearance of quiescent objects in our current sample.<fnref n="31" /> Source

#185 in PRIMER ($z = 4.51^{+0.16}_{-0.26}$, $\log(M_*/M_\odot) = 10.9$) is also picked as among the most reliable quiescent candidates by the NUVUVJ criterion ($P_{Q,50\%} = 0.79$). An entry with $z \sim 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. (2023b). However, we retrieve only #101962 (our #2876) above $z=4$ (Appendix B.4; see also Kocevski et al. 2023). Direct spectroscopic observations with JWST are necessary to break the current ceiling at $z \sim 4$ imposed by atmospheric hindering to ground-based telescopes and confirm the exact redshifts of these candidates.

## 4.6. Revisiting the Comparison with Simulations

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 IllustrisTNG 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 \ge 10^{10.6} M_\odot$ 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 QGs at $z>4$ would exacerbate the tension not only with these simulations but also with 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_\odot$ objects with sSFR $< 10^{-10} \text{ yr}^{-1}$ in EAGLE or IllustrisTNG at $z=4$ and above, while FLARES produces 2 dex lower number densities at $z=5$ ($n = 7.2 \times 10^{-8} Mpc^{-3}$; Lovell et al. 2022).

## 5. 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 UVJ selection and apply a novel technique based on Gaussian modeling of multiple colors—including an NUV 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.

1. We estimate $n \sim 2.5 \times 10^{-5} Mpc^{-3}$ for massive candidates ($\ge 10^{10.6} M_\odot$) with both selections, but substantial field-to-field variations of the order of $2-3 \times$. This is likely due to cosmic variance (~30%-50% uncertainty in the largest contiguous fields of ~30 arcmin² 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

<footnote n="31">³¹ Two bluer and potentially quenching objects are picked by the loosest UVJ selection at $z>5$. Their SEDs and color images are part of the overall release.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Robust $z > 4.5$ quiescent candidates. Left column: SEDs. Black squares and blue filled circles indicate the observed and best-fit photometry of each source. Light-gray squares mark observed flux densities with signal-to-noise ratio <3. Blue solid lines and shaded areas show the best-fit EAZY-PY models and their uncertainties, respectively. Middle 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'' \times 5''$.</figcaption></figure>

or quenching galaxies with consistent redshifts in the field with the highest number density.

2. We compile and homogenize the results of similar attempts to quantify the number densities of massive QGs 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 \sim (1-2) \times 10^{-5}$

$Mpc^{-3}$—not exceedingly far from what was established via ground-based observations.

3. 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 \sim 5$ are found in our systematic search and presented here.

4. 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} \le M_* < 10^{10.6} M_{\odot}$ similar to those at $\ge 10^{10.6} M_{\odot}$, 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 catalogs, images, and SED best-fit parameters) are publicly available. The continuous flow of new JWST imaging data and—soon—systematic spectroscopic coverage of large portions of the sky (e.g., Cosmos-Web, Casey et al. 2022; UNCOVER, Bezanson et al. 2022; GO 2665, PI: K. Glazebrook: Nanayakkara et al. 2022; 2362: PI: C. Marsan; 2285, PI: A. Carnall) will allow us to shrink the uncertainties due to cosmic variance and pursue the research avenues highlighted throughout the manuscript, starting with the necessary spectroscopic confirmation.

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 doi: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 Nos. 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 iMAGE (grant ANR-22-CE31-0007).

Facilities: JWST, HST.

Software: GRIZLI (Brammer & Matharu 2021; Brammer et al. 2022), EAZY-PY (Brammer et al. 2008), SEP (Barbary 2016), ASTROPY (Astropy Collaboration et al. 2022), ASTRODRIZZLE (Fruchter & Hook 2002; Koekemoer et al. 2003), GLAFIC (Oguri 2010, 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. Figure A1 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 A2.

#### A.1. CEERS

The Cosmic Evolution Early Release Science Survey (CEERS) is among the Director Discretionary Early Release Science (DD-ERS) programs (ERS 1345, PI: S. Finkelstein). It targeted the Extended Groth Strip (EGS) HST legacy field with several JWST instruments for imaging and, in the future, spectroscopy (Bagley et al. 2023 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).

<table><caption>**Table 3**
Filter Coverage in Each Field</caption><thead><tr><th>Field</th><th>JWST Wavelength<br>(μm)</th><th>HST Wavelength<br>(μm)</th></tr></thead><tbody><tr><td>CEERS</td><td>1.15, 1.5, 2, 2.7, 3.5, 4.1<sup>a</sup>, 4.4</td><td>0.6, 0.8, 0.44, 1.05, 1.25, 1.4, 1.6</td></tr><tr><td>Stephan's Quintet</td><td>1.5, 2, 2.7, 3.5, 4.4</td><td>...</td></tr><tr><td>PRIMER</td><td>0.9, 1.15, 1.5, 2, 2.7, 3.5, 4.1<sup>a</sup>, 4.4</td><td>0.44, 0.6, 0.8, 1.05, 1.25, 1.4, 1.6</td></tr><tr><td>NEP</td><td>0.9, 1.15, 1.5, 2, 2.7, 3.5, 4.1<sup>a</sup>, 4.4</td><td>0.44, 0.6</td></tr><tr><td>J1235</td><td>0.7, 0.9, 1.15, 1.5, 2, 2.7, 3.0<sup>a</sup>, 3.5, 4.4, 4.8<sup>a</sup></td><td>...</td></tr><tr><td>GLASS</td><td>0.9, 1.15, 1.5, 2, 2.7, 3.5, 4.4</td><td>0.44, 0.48, 0.6, 0.78, 0.8, 1.05, 1.25, 1.4, 1.6</td></tr><tr><td>Sunrise</td><td>0.9, 1.15, 1.5, 2, 2.7, 3.5, 4.1<sup>a</sup>, 4.4</td><td>0.44, 0.48, 0.6, 0.8, 1.05, 1.10, 1.25, 1.4, 1.6</td></tr><tr><td>SMACS 0723</td><td>0.9, 1.5, 2, 2.7, 3.5, 4.4</td><td>0.44, 0.6, 0.8, 1.05, 1.25, 1.4, 1.6</td></tr><tr><td>SGAS 1723</td><td>1.15, 1.5, 2, 2.7, 3.5, 4.4</td><td>0.48, 0.6, 0.78, 0.8, 1.05, 1.10, 1.4, 1.6</td></tr><tr><td>SPT 0418</td><td>1.15, 1.5, 2, 2.7, 3.5, 4.4</td><td>...</td></tr><tr><td>SPT 2147</td><td>2, 2.7, 3.5, 4.4</td><td>1.4</td></tr></tbody></table>

**Notes.** JWST NIRCam filter identifiers: Wide (W): 0.7 = F070W; 0.9 = F090W; 1.15 = F115W; 1.5 = F150W; 2 = F200W; 2.77 = F277W; 3.5 = F356W; 4.4 = F444W; Medium (M): 3.0 = F300M; 4.1 = F410M; 4.8 = F480M. HST filter identifiers: 0.44 = ACS/F435W; 0.48 = ACS/F475W; 0.6 = ACS/F606W; 0.78 = ACS/F775W; 0.8 = ACS/F814W; 1.05 = WFC3/F105W; 1.25 = WFC3/F125W; 1.4 = WFC3/F140W; 1.6 = WFC3/F160W.

<footnote><sup>a</sup> Medium-band filter.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure A1.** Observed NIRCam F444W magnitudes as a function of the photometric redshifts. Gray points indicate sources in each field as labeled. The color intensity scales as the density of points. Red circles show our *UVJ*-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. The complete figure set (5 elements) is available in the online Journal.</figcaption></figure>

(The complete figure set (5 images) is available.)
<figure><img src="image_11.png" /><figcaption>**Figure A2.** 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).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure A2. (Continued.)</figcaption></figure>
## 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 A2, 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².

## 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 MIRI. Here we considered the area covered with NIRCam in the UDS field, critically overlapping with the HST deep imaging from the Cosmic Assembly Near-infrared 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). The first spoke in the TDF was immediately released to the public. Here we considered the portion of the sky observed by NIRCam. Coverage with HST ACS/F435W and F606W is available.

## A.5. J1235

J1235 is a low ecliptic latitude field observed during commissioning with the largest compilation of wide and medium NIRCam filters in our collection in Cycle 0 (COM/NIRCam 1063, PI: B. Sunnquist). The goal was to verify to a 1% accuracy the flat-fielding after launch and to accumulate sky flats for future calibration programs. No HST imaging was available.

## A.6. GLASS Parallel

Parallel NIRCam observations were acquired while observing A2744 as part of the DD-ERS program "GLASS-JWST" (ERS 1324, PI: T. Treu; Treu et al. 2022). The parallel fields are sufficiently far from the cluster that gravitational lensing does not appreciably affect our work. A2744 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 WHL 0137-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. 2023)—even hosting a highly magnified star (Welch et al. 2022a). Being part of RELICS, ample HST ancillary data are 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, 2021). For an alternative modeling of the photometry on subgalactic spatial scales in this field and potential quiescent candidates at high redshift, we refer the reader to Abdurro'uf et al. (2023).

## A.8. SMACS 0723

SMACS 0723 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. SGAS 1723, SPT 0418, and SPT 2147

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 four 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. SPT 2147 was not imaged with the F115W and F150W filters.

# Appendix 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 B1 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-diameter 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 versions 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($\sigma_{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($\sigma_{MAD}$ = 0.30) and 0.13(0.24) dex, respectively. All things considered, the offsets are fully
<figure><img src="image_13.png" /></figure>

**Figure B1.** 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 open circles mark UVJ-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-diameter 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><img src="image_14.png" /></figure>

**Figure B2.** Effect of HST filter coverage. From left to right: rest-frame NUV, *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_{\text{phot}}$ in the top row and at fixed $z_{\text{phot}}$ (full filter set) in the bottom one to remove the effect of different redshifts. Gray circles indicate the whole galaxy sample. Blue circles mark bright objects with $F200W < 28$ mag in the redshift interval $0.3 \le z \le 6.5$ of interest, as labeled. Large filled circles indicate our selected QGs at $z > 3$ using the NUV/UVJ diagram with their uncertainties, color-coded according to the nominal $P_Q$ in the run with the full filter set.
<figure><img src="image_15.png" /><figcaption>**Figure B3.** 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. (2023b) (X-axes). Open and red filled circles indicate the full and “robust” samples of candidates identified by the authors in the CEERS field, respectively. We did not apply any correction to homogenize our Chabrier (2003) IMF with that of Kroupa (2001).</figcaption></figure>

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 UVJ- and NUVUVJ-selected QGs 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 NUV at $z > 3$. In Figure B2 we show the rest-frame NUV, U, V, and J flux densities in the two fitting runs, also removing the effect 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. Submillimetric Coverage and Spectroscopically Confirmed Objects

We cross-checked our list of candidate quiescent objects with available catalogs of submillimetric surveys in CEERS (450 and 850 µm down to $\sigma_{450} = 1.2$ and $\sigma_{850} = 0.2$ mJy beam$^{-1}$ with SCUBA-2 in the deep tier of the S2CLS survey, Zavala et al. 2017; $\sigma_{850} = 1.2$ mJy beam$^{-1}$ over the full survey, Geach et al. 2017), PRIMER (870 µm with ALMA from the AS2UDS survey targeting SCUBA-2 submillimeter galaxies from the S2CLS survey, Geach et al. 2017; and detecting sources as faint as 0.6 mJy at $>4.3\sigma$, Dudzevičiūtė et al. 2020; Cheng et al. 2023, based on a combination of archival data), and SMACS 0723 (1.1 mm observations with $\sigma_{1.1mm} = 66.1$ µJy beam$^{-1}$ from ALMA in the context of the ALCS Survey; Kokorev et al. 2022; Fujimoto et al. 2023). These limits correspond to SFR = 33–26 (CEERS/S2CLS-deep), 200–150 (S2CLS shallow); 25–19 (PRIMER/AS2UDS), 130–102 (S2CLS), and 25–16 $M_{\odot}$ yr$^{-1}$ (SMACS 0723/ALCS) at $z=3–6$, obtained by rescaling the $1\sigma$ rms with a modified blackbody with temperature $T_{dust} = 40$ K, $\beta=2$, $k_0 = 0.43$ cm$^2$ g$^{-1}$ at $\lambda_0 = 850$ µm (Li & Draine 2001), and accounting for the lesser effect of the cosmic microwave background (da Cunha et al. 2013). The S2CLS (shallow) survey covers all the CEERS and PRIMER fields. The deeper portion of the survey described in

Zavala et al. (2017) covers approximately 45% of our final samples in CEERS. The ALCS coverage of SMACS 0723 is of $\sim 3$ arcmin$^2$ centered on the cluster. AS2UDS and the ALMA archival observations are pointed and covered an area $\sim 600 \times$ 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 $\sim 5\sigma$ detection at 850 µm from SCUBA-2 at a 0.99 distance from a candidate UVJ QG 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 NUVUVJ criterion. However, several other possible optical/near-IR counterparts fall within the SCUBA-2 beam (S. Gillman et al. 2023, 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 both in our selections from Schreiber et al. (2018a) 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: ($\sim 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 (Valentino et al. 2022; G. Brammer et al. 2023, in preparation), nor in the 3D-HST survey (Skelton et al. 2014; Momcheva et al. 2016).

## B.4. Comparison with JWST-selected Photometric Quiescent Candidates in the Literature

Figure B3 shows the comparison between our F200W magnitudes and SED modeling results with EAZY-PY and those from Carnall et al. (2023b) 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 UVJ selection and that in Carnall et al. (2023b), especially for their “robust” sample. Sources #9844, 4921 in our catalog (78374, 76507 in Carnall et al. 2023b) are excluded by virtue of their blue colors, while #9131 (92564) has a large uncertainty on $V-J$ ($\sigma_{V-J} = 0.96$ mag). The overlap is less extended when imposing $P_{Q,50\%} \ge 0.1$. Sources below this threshold are at either the bluest (#9844, 4921) or the reddest end of the color distribution (e.g., #7432, 8556 = 40015, 42128), the latter being mainly occupied by dusty SFGs. We
note the fact that our photometry is extracted in 0.′′5 apertures and, thus, traces the properties of the central regions of galaxies. In the 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 B3). If any, our $z_{phot}$ seem to be systematically lower and $M_*$ larger than those derived by Carnall et al. (2023b) (Kocevski et al. 2023 also report lower redshift estimates). However, these offsets are in the realm of typical statistical and systematic uncertainties that different codes run with a variety of parameters can produce.

In addition, our selections do not retrieve the dusty candidate QGs at $z_{phot} \sim 5.4$ in SMACS 0723 presented in Rodighiero et al. (2023, ID#2 = KLAMA, #1536 in our catalog; R.A. = 110.70257564, decl. = -73.48472291 in the Gaia

DR3 astrometric reference). Our photometry and SED modeling place this object at $z_{phot} = 3.58_{-0.24}^{+0.60}$ and assign it values of $M_* = 3.0_{-0.8}^{+1.4} \times 10^{10} M_{\odot}$ and $P_{Q,50\%} \ll 0.1$.

We highlight the fact that the comparison with both these works is partially affected by the different JWST zero-point photometric calibration, an element in constant evolution to date.

## Appendix C: Stellar Mass Limits

Different mass limits could be a concern to draw comparison among fields with uneven photometric coverage and depth. In Figure C1, we show that our comparison 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

<figure><img src="image_16.png" /><figcaption>**Figure C1.** Stellar masses as a function of the photometric redshifts. Gray points indicate sources in each field as labeled. The color intensity scales as the density of points. The solid color lines mark the 90th 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 UVJ-selected red galaxies in the catalog at any redshifts (i.e., 90% of the UVJ-selected QGs at any redshifts lie above these lines in each bin). We show the 90th percentile of the distribution of galaxies in the CEERS catalog in each panel (dashed blue line). A direct comparison of the 90th percentiles is shown in the bottom right panel. For reference, the red circles show the location of our visually inspected UVJ-selected sample of quiescent candidates at $3 < z < 6.5$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure C2.** Stellar mass distributions of the selected samples. The solid lines mark the probability density functions of $M_*$ for the UVJ and NUVUVJ 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\%} \ge 0.7$ at $z > 4$. The red arrow marks its stellar mass estimate in the bottom right panel.</figcaption></figure>

adopted in Section 4 are well above the threshold where completeness is expected to drop. For reference, we show our visually inspected UVJ sample, more massive than the limit set by the 90th percentile of the $M_*$ distribution in the $3 < z < 6.5$ redshift interval in every field. This holds also for the robust NUVUVJ-selected sample.

Figure C2 shows the $M_*$ distributions of the UVJ and NUV–U, V–J selected sample. As noted in Section 3.3, candidates with $P_{Q,50\%} \ge 0.7$ maximally overlap with traditionally selected UVJ 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 UVJ pool of candidates do not introduce any immediately evident bias in mass.

## Appendix D: Literature Compilation

Table 4 includes complementary information on the comoving number densities of massive QGs at $3 \lesssim z \lesssim 4$ reported in Figure 5. The estimates are homogenized in terms of redshift interval and above the same mass limit of $\log(M_*/M_\odot) \gtrsim 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. (2018a) and Cecchi et al. (2019), while the redshift interval considered in Straatman et al. (2014) is $3.4 \lesssim 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. Finally, in Table 5 we report the comoving number densities of
<table><caption>**Table 4**
Comoving Number Densities of Quiescent Galaxies at $3 < z < 4$ in the Literature</caption><thead><tr><th>Label</th><th>n (Mpc<sup>-3</sup>)</th><th>Area/Box Side</th><th>Primary Selection</th><th>log(M<sub>☉</sub>) Limit (M<sub>☉</sub>)</th><th>IMF</th><th>References</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">Observations</td></tr><tr><td>Muzzin+13</td><td>2.7 × 10<sup>-6</sup></td><td>1.62 deg<sup>2</sup></td><td>UVJ (SMF)</td><td>10.6</td><td>Kroupa</td><td>M13, V20</td></tr><tr><td>Straatman+14</td><td>1.8<sup>+0.7</sup><sub>-0.7</sub> × 10<sup>-5</sup></td><td>121 arcmin<sup>2</sup></td><td>UVJ</td><td>10.6</td><td>Chabrier</td><td>S14</td></tr><tr><td>Davidzon+17</td><td>1.5<sup>+0.6</sup><sub>-0.7</sub> × 10<sup>-6</sup></td><td>0.62 deg<sup>2</sup></td><td>NUVrJ (SMF)</td><td>10.6</td><td>Chabrier</td><td>D17, V20</td></tr><tr><td rowspan="2">Schreiber+18</td><td>1.4<sup>+0.3</sup><sub>-0.3</sub> × 10<sup>-5</sup></td><td rowspan="2">442 arcmin<sup>2</sup></td><td>UVJ</td><td rowspan="2">10.5</td><td rowspan="2">Chabrier</td><td rowspan="2">S18</td></tr><tr><td>2.0<sup>+0.3</sup><sub>-0.3</sub> × 10<sup>-5</sup></td><td>sSFR</td></tr><tr><td rowspan="2">Merlin+19</td><td>4.9<sup>+0.6</sup><sub>-0.6</sub> × 10<sup>-6</sup></td><td rowspan="2">970 arcmin<sup>2</sup></td><td>SED (w/lines)</td><td rowspan="2">10.6</td><td rowspan="2">Chabrier</td><td rowspan="2">M18, M19</td></tr><tr><td>2.0<sup>+0.1</sup><sub>-0.1</sub> × 10<sup>-5</sup></td><td>SED (complete)</td></tr><tr><td>Cecchi+19</td><td>1.0<sup>+0.3</sup><sub>-0.3</sub> × 10<sup>-6</sup></td><td>1.38 deg<sup>2</sup></td><td>NUVrJ+sSFR</td><td>10.5</td><td>Chabrier</td><td>C19, This work</td></tr><tr><td>Girelli+19</td><td>1.5<sup>+0.5</sup><sub>-0.7</sub> × 10<sup>-6</sup></td><td>1.38 deg<sup>2</sup></td><td>Observed colors</td><td>10.6</td><td>Chabrier</td><td>G19, V20</td></tr><tr><td>Shahidi+20</td><td>4.9<sup>+1.6</sup><sub>-1.2</sub> × 10<sup>-6</sup></td><td>964 arcmin<sup>2</sup></td><td>Balmer, UVJ, SED</td><td>10.6</td><td>Chabrier</td><td>S20, This work</td></tr><tr><td rowspan="2">Carnall+20</td><td>1.1<sup>+0.4</sup><sub>-0.3</sub> × 10<sup>-5</sup></td><td rowspan="2">370 arcmin<sup>2</sup></td><td>sSFR (Full)</td><td rowspan="2">10.6</td><td rowspan="2">Kroupa</td><td rowspan="2">C20, This work</td></tr><tr><td>5.1<sup>+3.0</sup><sub>-2.0</sub> × 10<sup>-6</sup></td><td>sSFR (Robust)</td></tr><tr><td>Weaver+22</td><td>9.4<sup>+1.7</sup><sub>-0.7</sub> × 10<sup>-6</sup></td><td>1.27 deg<sup>2</sup></td><td>NUVrJ</td><td>10.6</td><td>Chabrier</td><td>W22, This work</td></tr><tr><td>Gould+23</td><td>1.2<sup>+0.4</sup><sub>-0.4</sub> × 10<sup>-5</sup></td><td>1.27 deg<sup>2</sup></td><td>GMM</td><td>10.6</td><td>Chabrier</td><td>G23</td></tr><tr><td rowspan="2">Carnall+23</td><td>6.3<sup>+3.8</sup><sub>-2.5</sub> × 10<sup>-5</sup></td><td rowspan="2">30 arcmin<sup>2</sup></td><td>sSFR (Full)</td><td rowspan="2">10.6</td><td rowspan="2">Kroupa</td><td rowspan="2">C23, This work</td></tr><tr><td>4.2<sup>+3.3</sup><sub>-2.0</sub> × 10<sup>-5</sup></td><td>sSFR (Robust)</td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">Simulations</td></tr><tr><td rowspan="2">Illustris-1</td><td>3.3<sup>+0.8</sup><sub>-0.8</sub> × 10<sup>-6</sup></td><td rowspan="2">107 cMpc</td><td>sSFR (z = 3.0)</td><td rowspan="2">10.6</td><td rowspan="2">Chabrier</td><td rowspan="2">V20</td></tr><tr><td>&lt;8.1 × 10<sup>-7</sup></td><td>sSFR (z = 3.7)</td></tr><tr><td rowspan="2">IllustrisTNG100</td><td>5.4<sup>+0.0</sup><sub>-1.0</sub> × 10<sup>-5</sup></td><td rowspan="2">107 cMpc</td><td>sSFR (z = 3.0)</td><td rowspan="2">10.6</td><td rowspan="2">Chabrier</td><td rowspan="2">V20</td></tr><tr><td>7.8 × 10<sup>-6</sup></td><td>sSFR (z = 3.7)</td></tr><tr><td rowspan="2">IllustrisTNG300</td><td>3.0<sup>+0.0</sup><sub>-0.5</sub> × 10<sup>-5</sup></td><td rowspan="2">293 cMpc</td><td>sSFR (z = 3.0)</td><td rowspan="2">10.6</td><td rowspan="2">Chabrier</td><td rowspan="2">V20</td></tr><tr><td>2.6<sup>+0.0</sup><sub>-0.5</sub> × 10<sup>-6</sup></td><td>sSFR (z = 3.7)</td></tr><tr><td rowspan="2">EAGLE</td><td>1.0 × 10<sup>-6</sup></td><td rowspan="2">100 cMpc</td><td>sSFR (z = 3.0)</td><td rowspan="2">10.6</td><td rowspan="2">Chabrier</td><td rowspan="2">This work</td></tr><tr><td>&lt;1.8 × 10<sup>-6</sup></td><td>sSFR (z = 3.9)</td></tr></tbody></table>

**Note.** Primary selection: UVJ, NUVrJ = rest-frame color diagrams; SMF = integration of stellar mass functions (from best-fit parameters); sSFR = cuts in sSFRs estimated via SED modeling; SED = composite cuts based on SED modeling; observed colors; Balmer = Balmer break; GMM = Gaussian mixture modeling in the NUV - U, U - V, V - J space. IMF: Kroupa (2001), Chabrier (2003). We do not apply any mass correction to convert one IMF to the other.

**References.** M13 = Muzzin et al. (2013); S14 = Straatman et al. (2014); D17 = Davidzon et al. (2017); S18 = Schreiber et al. (2018a); M18 = Merlin et al. (2018); M19 = Merlin et al. (2019); C19 = Cecchi et al. (2019); G19 = Girelli et al. (2019); V20 = Valentino et al. (2020); S20 = Shahidi et al. (2020); C20 = Carnall et al. (2020); W22 = Weaver et al. (2022a); G23 = Gould et al. (2023); C23 = Carnall et al. (2023b).
<table><caption>**Table 5**
Comoving Number Densities of Quiescent Galaxies at $3 < z < 4$</caption><thead><tr><th>Field</th><th>Area<br>(arcmin<sup>2</sup>)</th><th>UVJ<br>Strict</th><th>UVJ<br>Padded</th><th>NUVUVJ<br>P<sub>Q,50%</sub></th><th>σ<sub>CV</sub><br>(%)</th></tr></thead><tbody><tr><td>CEERS</td><td>34.7</td><td>3.2<sup>+2.8</sup><sub>-1.6</sub><br>6.6<sup>+3.5</sup><sub>-2.4</sub></td><td>4.0<sup>+3.0</sup><sub>-1.8</sub><br>7.5<sup>+3.6</sup><sub>-2.6</sub></td><td>3.9<sup>+2.9</sup><sub>-1.8</sub><br>7.6<sup>+3.6</sup><sub>-2.6</sub></td><td>0.29<br>0.53</td></tr><tr><td>Stephan's Quintet</td><td>35.0</td><td>6.4<sup>+3.4</sup><sub>-2.4</sub><br>0.0<sup>+1.7</sup><sub>-0.0</sub></td><td>7.4<sup>+3.6</sup><sub>-2.5</sub><br>0.0<sup>+1.8</sup><sub>-0.0</sub></td><td>3.9<sup>+2.9</sup><sub>-1.8</sub><br>0.0<sup>+1.7</sup><sub>-0.0</sub></td><td>0.30<br>0.55</td></tr><tr><td>PRIMER</td><td>21.9</td><td>2.9<sup>+3.8</sup><sub>-1.9</sub><br>3.4<sup>+4.0</sup><sub>-2.1</sub></td><td>3.0<sup>+3.8</sup><sub>-1.9</sub><br>4.5<sup>+4.3</sup><sub>-2.4</sub></td><td>1.6<sup>+3.4</sup><sub>-1.3</sub><br>3.4<sup>+4.0</sup><sub>-2.1</sub></td><td>0.32<br>0.58</td></tr><tr><td>NEP</td><td>9.7</td><td>7.5<sup>+8.9</sup><sub>-4.6</sub><br>0.0<sup>+6.0</sup><sub>-0.0</sub></td><td>7.5<sup>+8.9</sup><sub>-4.6</sub><br>0.0<sup>+6.6</sup><sub>-0.0</sub></td><td>5.0<sup>+8.1</sup><sub>-3.6</sub><br>0.0<sup>+6.0</sup><sub>-0.0</sub></td><td>0.34<br>0.61</td></tr><tr><td>J1235</td><td>9.0</td><td>0.3<sup>+6.7</sup><sub>-0.3</sub><br>0.0<sup>+6.5</sup><sub>-0.0</sub></td><td>0.3<sup>+6.7</sup><sub>-0.3</sub><br>0.3<sup>+6.6</sup><sub>-0.3</sub></td><td>3.7<sup>+8.2</sup><sub>-3.0</sub><br>0.3<sup>+6.6</sup><sub>-0.3</sub></td><td>0.34<br>0.62</td></tr><tr><td>GLASS</td><td>8.5</td><td>0.0<sup>+6.9</sup><sub>-0.0</sub><br>0.0<sup>+6.9</sup><sub>-0.0</sub></td><td>0.0<sup>+6.9</sup><sub>-0.0</sub><br>0.0<sup>+6.9</sup><sub>-0.0</sub></td><td>0.0<sup>+6.9</sup><sub>-0.0</sub><br>0.0<sup>+6.9</sup><sub>-0.0</sub></td><td>0.34<br>0.62</td></tr><tr><td>Sunrise</td><td>7.3</td><td>0.0<sup>+8.1</sup><sub>-0.0</sub><br>0.0<sup>+8.1</sup><sub>-0.0</sub></td><td>0.0<sup>+8.1</sup><sub>-0.0</sub><br>0.0<sup>+8.1</sup><sub>-0.0</sub></td><td>0.0<sup>+8.1</sup><sub>-0.0</sub><br>0.0<sup>+8.1</sup><sub>-0.0</sub></td><td>0.34<br>0.62</td></tr><tr><td>SMACS 0723</td><td>6.5</td><td>6.2<sup>+11.7</sup><sub>-4.7</sub><br>0.8<sup>+9.5</sup><sub>-0.8</sub></td><td>6.2<sup>+11.7</sup><sub>-4.7</sub><br>0.8<sup>+9.5</sup><sub>-0.8</sub></td><td>4.0<sup>+10.9</sup><sub>-3.5</sub><br>0.0<sup>+9.0</sup><sub>-0.0</sub></td><td>0.34<br>0.63</td></tr><tr><td>SGAS 1723</td><td>5.3</td><td>5.4<sup>+13.5</sup><sub>-4.6</sub><br>6.1<sup>+13.8</sup><sub>-5.0</sub></td><td>1.9<sup>+12.0</sup><sub>-1.9</sub><br>6.1<sup>+13.8</sup><sub>-5.0</sub></td><td>0.0<sup>+11.0</sup><sub>-0.0</sub><br>0.0<sup>+11.0</sup><sub>-0.0</sub></td><td>0.35<br>0.64</td></tr><tr><td>SPT 0418</td><td>5.0</td><td>2.5<sup>+13.0</sup><sub>-2.5</sub><br>0.0<sup>+11.7</sup><sub>-0.0</sub></td><td>2.5<sup>+13.0</sup><sub>-2.5</sub><br>0.0<sup>+11.7</sup><sub>-0.0</sub></td><td>0.0<sup>+11.7</sup><sub>-0.0</sub><br>0.0<sup>+11.7</sup><sub>-0.0</sub></td><td>0.35<br>0.65</td></tr><tr><td>SPT 2147</td><td>2.3</td><td>0.0<sup>+25.4</sup><sub>-0.0</sub><br>0.0<sup>+25.4</sup><sub>-0.0</sub></td><td>0.0<sup>+25.4</sup><sub>-0.0</sub><br>0.0<sup>+25.4</sup><sub>-0.0</sub></td><td>0.0<sup>+25.4</sup><sub>-0.0</sub><br>0.0<sup>+25.4</sup><sub>-0.0</sub></td><td>0.37<br>0.67</td></tr><tr><td>Combined</td><td>145.1</td><td>3.9<sup>+1.2</sup><sub>-0.9</sub><br>2.4<sup>+1.0</sup><sub>-0.7</sub></td><td>4.1<sup>+1.2</sup><sub>-0.9</sub><br>2.7<sup>+1.0</sup><sub>-0.8</sub></td><td>2.8<sup>+1.0</sup><sub>-0.8</sub><br>2.3<sup>+1.0</sup><sub>-0.7</sub></td><td>0.10<br>0.18</td></tr></tbody></table>

Note. The comoving number densities are expressed in units of $10^{-5}$ Mpc$^{-3}$. Each field has two entries: the first and second rows refer to the log($M_*/M_c$) = [9.5, 10.6) and $\ge$10.6 bins, respectively. The uncertainties reflect the Poissonian $1\sigma$ confidence interval. Upper limits are at $1\sigma$ using the same approach (Gehrels 1986). Statistical uncertainties are accounted for by integrating the $p(z)$ within $3 < z < 4$. The uncertainties due to cosmic variance are expressed as fractional $\sigma_{CV}$ deviations (Section 4.1). The selections are described in Section 3. The adopted threshold for the NUVUVJ selection is $P_{Q,50\%} \ge 0.1$.

quiescent objects in the same redshift and mass bins as in Table 2, but for each individual field that we considered in this work.

## ORCID iDs

Francesco Valentino https://orcid.org/0000-0001-6477-4011

Gabriel Brammer https://orcid.org/0000-0003-2680-005X

Katriona M. L. Gould https://orcid.org/0000-0003-4196-5960

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

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

Christian Kragh Jespersen https://orcid.org/0000-0002-8896-6496

Aswin P. Vijayan https://orcid.org/0000-0002-1905-4194

John R. Weaver https://orcid.org/0000-0003-1614-196X

Kei Ito https://orcid.org/0000-0002-9453-0381

Masayuki Tanaka https://orcid.org/0000-0002-5011-5178

Olivier Ilbert https://orcid.org/0000-0002-7303-4397

Georgios E. Magdis https://orcid.org/0000-0002-4872-2294

Katherine E. Whitaker https://orcid.org/0000-0001-7160-3632

Andreas L. Faisst https://orcid.org/0000-0002-9382-9832

Anna Gallazzi https://orcid.org/0000-0002-9656-1800

Steven Gillman https://orcid.org/0000-0001-9885-4589

Clara Giménez-Arteaga https://orcid.org/0000-0001-9419-9505

Carlos Gómez-Guijarro https://orcid.org/0000-0002-4085-9165

Mariko Kubo https://orcid.org/0000-0002-7598-5292

Kasper E. Heintz https://orcid.org/0000-0002-9389-7413

Michaela Hirschmann https://orcid.org/0000-0002-3301-3321

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

Masato Onodera https://orcid.org/0000-0003-3228-7264

Francesca Rizzo https://orcid.org/0000-0001-9705-2461

Minju Lee https://orcid.org/0000-0002-2419-3068

Victoria Strait https://orcid.org/0000-0002-6338-7295

Sune Toft https://orcid.org/0000-0003-3631-7176

## References

Abdurro'uf, C. D., Jung, I., Ferguson, H. C., et al. 2023, ApJ, 945, 117

Akhshik, M., Whitaker, K. E., Leja, J., et al. 2023, ApJ, 943, 179

Antwi-Danso, J., Papovich, C., Leja, J., et al. 2023, ApJ, 943, 166

Arnouts, S., Le Floc'h, E., Chevallard, J., et al. 2013, A&A, 558, A67

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

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

Barbary, K. 2016, JOSS, 1, 58

Barrufet, L., Oesch, P. A., Weibel, A., et al. 2022, arXiv:2207.14733

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

Bezanson, R., Labbe, I., Whitaker, K. E., et al. 2022, arXiv:2212.04026

Boyer, M. L., Anderson, J., Gennaro, M., et al. 2022, RNAAS, 6, 191

Bradley, L. D., Coe, D., Brammer, G., et al. 2022, arXiv:2210.01777

Brammer, G. 2022, Preliminary updates to the NIRCam photometric calibration, v1, Zenodo, doi:10.5281/zenodo.7143382

Brammer, G. 2023, JWST image mosaics grizli-v4, University of Copenhagen, doi:10.17894/UCPH.E3D897AF-233A-4F01-A893-7B0FAD1F66C2

Brammer, G., & Matharu, J. 2021, gbrammer/grizli: Release 2021, v1.3.2, Zenodo, doi:10.5281/zenodo.1146904

Brammer, G., Strait, V., Matharu, J., & Momcheva, I. 2022, grizli, v1.5.0, Zenodo, doi:10.5281/zenodo.6672538

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

Brammer, G. B., Whitaker, K. E., van Dokkum, P. G., et al. 2011, ApJ, 739, 24

Byler, N. 2018, cloudyFSPS, v1.0.0, Zenodo, doi:10.5281/zenodo.1156412

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

Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019, ApJ, 873, 44

Carnall, A. C., McLeod, D. J., McLure, R. J., et al. 2023b, MNRAS, 520, 3974

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2023c, arXiv:2301.11413

Carnall, A. C., Walker, S., McLure, R. J., et al. 2020, MNRAS, 496, 695

Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2022, arXiv:2211.07865

Cecchi, R., Bolzonella, M., Cimatti, A., & Girelli, G. 2019, ApJL, 880, L14

Chabrier, G. 2003, PASP, 115, 763

Cheng, C., Huang, J.-S., Smail, I., et al. 2023, ApJL, 942, L19

Chuang, C.-Y., & Lin, Y.-T. 2023, ApJ, 944, 207

Ciesla, L., Boselli, A., Elbaz, D., et al. 2016, A&A, 585, A43

Cimatti, A., Cassata, P., Pozzetti, L., et al. 2008, A&A, 482, 21

Coe, D., Salmon, B., Bradač, M., et al. 2019, ApJ, 884, 85

Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

De Lucia, G., Springel, V., White, S. D. M., Croton, D., & Kauffmann, G. 2006, MNRAS, 366, 499

D'Eugenio, C., Daddi, E., Gobat, R., et al. 2020, ApJL, 892, L2

D'Eugenio, C., Daddi, E., Gobat, R., et al. 2021, A&A, 653, A32

Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485, 4817

Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828

Esdaile, J., Glazebrook, K., Labbé, I., et al. 2021, ApJL, 908, L35
Faisst, A. L., Chary, R. R., Brammer, G., & Toft, S. 2022, ApJL, 941, L11

Fontana, A., Santini, P., Grazian, A., et al. 2009, A&A, 501, 15

Forrest, B., Annunziatella, M., Wilson, G., et al. 2020a, ApJL, 890, L1

Forrest, B., Marsan, Z. C., Annunziatella, M., et al. 2020b, ApJ, 903, 47

Forrest, B., Wilson, G., Muzzin, A., et al. 2022, ApJ, 938, 109

Fruchter, A. S., & Hook, R. N. 2002, PASP, 114, 144

Fujimoto, S., Finkelstein, S. L., Burgarella, D., et al. 2022, arXiv:2211.03896

Fujimoto, S., Kohno, K., Ouchi, M., et al. 2023, arXiv:2303.01658

Fujimoto, S., Ouchi, M., Ono, Y., et al. 2016, ApJS, 222, 1

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1

Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789

Gehrels, N. 1986, ApJ, 303, 336

Giménez-Arteaga, C., Oesch, P. A., Brammer, G. B., et al. 2022, arXiv:2212.08670

Girelli, G., Bolzonella, M., & Cimatti, A. 2019, A&A, 632, A80

Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Natur, 544, 71

Gould, K. M. L., Brammer, G., Valentino, F., et al. 2023, arXiv:2302.10934

Green, G. 2018, JOSS, 3, 695

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5

Ito, K., Tanaka, M., Miyaji, T., et al. 2022, ApJ, 929, 53

Ito, K., Tanaka, M., Valentino, F., et al. 2023, ApJL, 945, L9

Iyer, K. G., Gawiser, E., Faber, S. M., et al. 2019, ApJ, 879, 116

Jespersen, C. K., Cranmer, M., Melchior, P., et al. 2022, ApJ, 941, 7

Kalita, B. S., Daddi, E., D'Eugenio, C., et al. 2021, ApJL, 917, L17

Kocevski, D. D., Barro, G., McGrath, E. J., et al. 2023, ApJL, 946, L14

Koekemoer, A. M., Fruchter, A. S., Hook, R. N., & Hack, W. 2003, in HST Calibration Workshop: Hubble after the Installation of the ACS and the NICMOS Cooling System, ed. S. Arribas, A. Koekemoer, & B. Whitmore (Baltimore, MD: STScI), 337

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

Kriek, M., & Conroy, C. 2013, ApJL, 775, L16

Kron, R. G. 1980, ApJS, 43, 305

Kroupa, P. 2001, MNRAS, 322, 231

Kubo, M., Tanaka, M., Yabe, K., et al. 2018, ApJ, 867, 1

Kubo, M., Umehata, H., Matsuda, Y., et al. 2021, ApJ, 919, 6

Kubo, M., Umehata, H., Matsuda, Y., et al. 2022, ApJ, 935, 89

Labbe, I., van Dokkum, P., Nelson, E., et al. 2022, arXiv:2207.12446

Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599

Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019a, ApJ, 876, 3

Leja, J., Tacchella, S., & Conroy, C. 2019b, ApJL, 880, L9

Li, A., & Draine, B. T. 2001, ApJ, 554, 778

Long, A. S., Casey, C. M., Lagos, C. d. P., et al. 2022, arXiv:2211.02072

Lovell, C. C., Roper, W., Vijayan, A. P., et al. 2022, arXiv:2211.07540

Lovell, C. C., Vijayan, A. P., Thomas, P. A., et al. 2021, MNRAS, 500, 2127

Lustig, P., Strazzullo, V., D'Eugenio, C., et al. 2021, MNRAS, 501, 2659

Maíz Apellániz, J. 2006, AJ, 131, 1184

Man, A., & Belli, S. 2018, NatAs, 2, 695

Manning, S. M., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 925, 23

Marchesini, D., Brammer, G., Morishita, T., et al. 2023, ApJL, 942, L25

Marsan, Z. C., Marchesini, D., Brammer, G. B., et al. 2015, ApJ, 801, 133

Marsan, Z. C., Marchesini, D., Brammer, G. B., et al. 2017, ApJ, 842, 21

Marsan, Z. C., Muzzin, A., Marchesini, D., et al. 2022, ApJ, 924, 25

Mawatari, K., Inoue, A. K., Hashimoto, T., et al. 2020, ApJ, 889, 137

McConachie, I., Wilson, G., Forrest, B., et al. 2022, ApJ, 926, 37

Mehta, V., Scarlata, C., Capak, P., et al. 2018, ApJS, 235, 36

Merlin, E., Fontana, A., Castellano, M., et al. 2018, MNRAS, 473, 2098

Merlin, E., Fortuni, F., Torelli, M., et al. 2019, MNRAS, 490, 3309

Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27

Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJS, 206, 8

Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJL, 699, L178

Nanayakkara, T., Glazebrook, K., Jacobs, C., et al. 2022, arXiv:2212.11638

Nardiello, D., Bedin, L. R., Burgasser, A., et al. 2022, MNRAS, 517, 484

Nelson, D., Springel, V., Pillepich, A., et al. 2019, ComAC, 6, 2

Nelson, E. J., Suess, K. A., Bezanson, R., et al. 2022, arXiv:2208.01630

Oguri, M. 2010, PASJ, 62, 1017

Oguri, M. 2021, PASP, 133, 074504

Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, JMLR, 12, 2825

Peng, Y., Maiolino, R., & Cochrane, R. 2015, Natur, 521, 192

Pérez-González, P. G., Barro, G., Annunziatella, M., et al. 2022, arXiv:2211.00045

Pontoppidan, K. M., Barrientes, J., Blome, C., et al. 2022, ApJL, 936, L14

Rodighiero, G., Bisigello, L., Iani, E., et al. 2023, MNRAS Lett., 518, L19

Salmon, B., Coe, D., Bradley, L., et al. 2020, ApJ, 889, 189

Santini, P., Castellano, M., Fontana, A., et al. 2022, ApJ, 940, 135

Santini, P., Castellano, M., Merlin, E., et al. 2021, A&A, 652, A30

Santini, P., Merlin, E., Fontana, A., et al. 2019, MNRAS, 486, 560

Saracco, P., Marchesini, D., La Barbera, F., et al. 2020, ApJ, 905, 40

Schawinski, K., Urry, C. M., Simmons, B. D., et al. 2014, MNRAS, 440, 889

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Schlaflay, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018a, A&A, 618, A85

Schreiber, C., Labbé, I., Glazebrook, K., et al. 2018b, A&A, 611, A22

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Shahidi, A., Mobasher, B., Nayyeri, H., et al. 2020, ApJ, 897, 44

Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Spitler, L. R., Straatman, C. M. S., Labbé, I., et al. 2014, ApJ, 787, L36

Stach, S. M., Dudzevičiūtė, U., Smail, I., et al. 2019, MNRAS, 487, 4648

Steinhardt, C. L., Jespersen, C. K., & Linzer, N. B. 2021, ApJ, 923, 8

Straatman, C. M. S., Labbé, I., Spitler, L. R., et al. 2014, ApJL, 783, L14

Strazzullo, V., Daddi, E., Gobat, R., et al. 2015, A&A, 576, L6

Suzuki, T. L., Glazebrook, K., Schreiber, C., et al. 2022, ApJ, 936, 61

Tanaka, A., Valentino, F., Toft, S., et al. 2019, ApJL, 885, L34

Toft, S., Smolčić, V., Magnelli, B., et al. 2014, ApJ, 782, 68

Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110

Valentino, F., Brammer, G., Fujimoto, S., et al. 2022, ApJL, 929, L9

Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93

Vanzella, E., Claeyssens, A., Welch, B., et al. 2023, ApJ, 945, 53

Vijayan, A. P., Lovell, C. C., Wilkins, S. M., et al. 2021, MNRAS, 501, 3289

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Natur, 572, 211

Weaver, J. R., Davidzon, I., Toft, S., et al. 2022a, arXiv:2212.02512

Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022b, ApJS, 258, 11

Welch, B., Coe, D., Diego, J. M., et al. 2022a, Natur, 603, 815

Welch, B., Coe, D., Zackrisson, E., et al. 2022b, ApJL, 940, L1

Williams, C. C., Labbe, I., Spilker, J., et al. 2019, ApJ, 884, 154

Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879

Zavala, J. A., Aretxaga, I., Geach, J. E., et al. 2017, MNRAS, 464, 3369