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

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 002 pixel scale for the Short Wavelength (SW) NIRCam bands and to 004 for all the remaining JWST and Hubble Space Telescope (HST) filters. We provide further details about the individual fields in Appendix A.

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 05 and correct to the "total" values within an elliptical Kron aperture (Kron 1980). ²³  The aperture correction is computed on the LW detection image and applied to all bands. The depths in the reference 05 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). ²⁴

We utilize eazy-py ²⁵  (Brammer et al. 2008) to estimate photometric redshifts, rest-frame colors, and stellar masses from the 05-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. ²⁶  The templates are created adopting a Chabrier (2003) IMF and applying a Kriek & Conroy (2013) dust attenuation law (dust index δ = − 0.1, R _V = 3.1), where the maximum allowed attenuation is also redshift dependent. A fixed grid of nebular emission lines and continuum from cloudy (v13.03) models is added to the templates within fsps (metallicity: $\mathrm{log}(Z/{Z}_{\odot })\in [-1.2,0]$, ionization parameter $\mathrm{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. ²⁷  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 σ_NMAD = 0.0268 (0.0187) for the spectroscopic sample (clipping catastrophic outliers). ²⁸ 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.

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 (λ = 2800 Å), 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.

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σ uncertainties on the colors consistent with the selection box as long as σ_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.

Zoom In Zoom Out Reset image size

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σ 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.

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 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 × 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%} ²⁹  and its 1σ uncertainties. We also list the nominal P _Q associated with the best-fit colors in our catalogs.

In the rest of our analysis, we adopt a cut at P _{Q ,50%} ≥ 0.1 to select candidate 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 ≤ P _{Q ,50%} < 0.7 and 18/20 (∼90%) truly passive galaxy candidates with P _{Q ,50%} ≥ 0.7. Their location in the projected NUV − U, V − J plane is shown in Figure 2.

Zoom In Zoom Out Reset image size

As noted in Gould et al. (2023), a selection in the NUVUVJ arguably outperforms the classical UVJ in selecting quiescent (passive and recently quenched or post-starbust) galaxies at z > 3. However, the two criteria partially overlap and identify the same quiescent sources—to an extent fixed by the P _Q threshold and the exact location of the selection box in the 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 NUVUVJ 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%} ≥ 0.7 and 0.1 ≤ P _{Q ,50%} < 0.7 are part of the UVJ sample. Moreover, 16/18 objects with P _{Q ,50%} ≥ 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.

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

3.4.2. Availability of HST Imaging

3.4.3. Dusty Star-forming or High-redshift Contaminants

3.4.4. Spectroscopically Confirmed and Alternative JWST-based Photometric Samples

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 $\mathrm{log}({M}_{\star }/{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 cosmic variance for the bin at $\mathrm{log}({M}_{\star }/{M}_{\odot })=9.5\mbox{--}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

$$\begin{eqnarray}&&{\sigma }_{\mathrm{CV},\mathrm{total}}=\sqrt{\displaystyle \frac{1}{\displaystyle \sum _{\mathrm{fields}}{\sigma }_{\mathrm{CV},\mathrm{field}}^{-2}}},\end{eqnarray} \tag{ 1 }$$

which assumes that all fields are independent. ³⁰  The relative uncertainties due to cosmic variance in the combined field are reported in Table 2.

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 simulated galaxies with sSFR ≤ 10⁻¹⁰ yr⁻¹ 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).

Zoom In Zoom Out Reset image size

Our number density estimates from the combined fields are of the order of ∼2.5 × 10⁻⁵ Mpc⁻³, consistent with some of the determinations with similar color or sSFR cuts (Schreiber et al. 2018a; Merlin et al. 2019). Our estimates are ∼2× larger than the most recent measurements in the largest contiguous survey among those considered, COSMOS (Weaver et al. 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.

We notice a substantial field-to-field variation especially when focusing on the most massive galaxies. Compared with the average number density on the full combined field of ∼145 arcmin², we find per-field value oscillations of a factor of 2–3× (Table 5). We ascribe these differences to cosmic variance and to the fact that massive quiescent systems might already be signposts of distant overdensities and protoclusters—massive halos able to fast-track galaxy evolution. In Table 5, we report the fractional 1σ uncertainties due to cosmic variance in each field. Taken individually, an uncertainty of ∼30%–50% affects the number densities in the two mass bins for the largest contiguous areas that we considered. The Stephan's Quintet and CEERS fields are emblematic in this sense, appearing under- and overdense despite a similar sky coverage. CEERS displays the largest number densities of QGs >10^10.6 M_⊙ 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 ∼ 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 ≳ 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.

Figures 3 and 4 show the number densities at lower stellar masses (10^9.5 M_⊙ ≤ M_⋆ < 10^10.6 M_⊙), now safely accessible with JWST also at such high redshifts. In fact, the low-mass end of the M_⋆ distributions starts declining at thresholds as low as ∼10⁸ M_⊙ at 3 < z < 6.5 (excluding SPT 2147; see Figure C1). The lower limit of M_⋆ = 10^9.5 M_⊙ chosen for the calculation is similar to that fixed in some of the works listed in Table 4. Thus, it allows us to derive a relatively straightforward comparison, discounting some of the systematics mentioned above. At 3 < z < 4, we estimate modestly higher (∼1.2–1.6 ×) number densities than in the most massive bin, but in agreement within the uncertainties. The difference with previous works is similar in each bin, when available. This is consistent with the expected shape of the stellar mass function of red QGs, roughly peaking and flattening or turning over at ∼10^10.6 M_⊙ 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 ∼ 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.

We can now look at galaxies at 4 < z < 5. As in the 3 < z < 4 interval, when considering the combined fields, we estimate number densities that are ∼2.5–4.5 × smaller above and below M_⋆ = 10^10.6 M_⊙ 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 ∼2.5–1 × when we consider only the same field and the "robust" sample in their work. This seems to suggest that cosmic variance and early overdense environments are effective in producing substantial field-to-field variations also at z > 4 (Section 4.3). For reference, our number density estimates in the same massive bin (M_⋆ ≥ 10^10.6 M_⊙) are consistent with those in the COSMOS field by Weaver et al. (2022a), but 1.8× 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_⊙, 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 ± 0.3) × 10⁻⁵ Mpc⁻³ for M_⋆ > 10^9.9 M_⊙; Weaver et al. 2022a) and 1.6–2× larger than in large-field HST surveys such as CANDELS (∼7.9 × 10⁻⁶ Mpc⁻³ for the "complete" sample at M_⋆ > 5 × 10⁹ M_⊙; Merlin et al. 2019). Finally, the upper limits on number densities for the highest-redshift bin at 5 < z < 6.5 should be taken with caution, given the area covered in our analysis (Table 2).

Focusing on the highest envelope of the redshift interval spanned by our sample, we find a few promising candidates at z > 4.5. The SEDs and three-color images of the 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 ≳ 5.2, which signposts the earliest epoch of appearance of quiescent objects in our current sample. ³¹  Source #185 in PRIMER ($z={4.51}_{-0.26}^{+0.16}$, $\mathrm{log}({M}_{\star }/{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 ∼ 3.2 at <03 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}_{\mathrm{phot}}={2.78}_{-0.08}^{+0.10}$). 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 ∼ 4 imposed by atmospheric hindering to ground-based telescopes and confirm the exact redshifts of these candidates.

Zoom In Zoom Out Reset image size

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 ≥ 10^10.6 M_⊙ QGs already at z = 3, while the situation seems partially alleviated if one includes lower-mass galaxies in the calculation (Merlin et al. 2019; Lovell et al. 2022). Spectroscopically confirmed massive 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_⊙ objects with sSFR < 10⁻¹⁰ yr⁻¹ in EAGLE or IllustrisTNG at z = 4 and above, while FLARES produces 2 dex lower number densities at z = 5 (n = 7.2 × 10⁻⁸ Mpc⁻³; Lovell et al. 2022).

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).

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².

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).

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.

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.

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.

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).

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.

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.

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 <05 separation. Figure B1 shows the comparison in z_phot, M_⋆, total HST/F160W photometry computed from our reference 05-diameter aperture, and that within a common 07-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 05-diameter apertures considered here are fainter than those from 07 in 3D-HST in CEERS and PRIMER: the median differences are 0.11(σ_MAD = 0.14) and 0.16(0.14) mag, respectively. However, at fixed aperture, the total magnitudes are fully consistent. The difference arises from the detection bands and where the aperture correction is computed: F160W for 3D-HST and the NIRCam combined LW image in our analysis. Finally, the total M_⋆ are systematically lower in 3D-HST than in our JWST-based catalogs of CEERS and PRIMER, with median differences and MAD of 0.19(σ_MAD = 0.30) and 0.13(0.24) dex, respectively. All things considered, the offsets are fully ascribable to the different recipes adopted to estimate these quantities and consistent with typical systematic uncertainties inevitably present when we compare different catalogs of the same sources. Our samples of UVJ- and NUVUVJ-selected QGs at z > 3 do not appreciably deviate from these trends.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

We cross-checked our list of candidate quiescent objects with available catalogs of submillimetric surveys in CEERS (450 and 850 μm down to σ₄₅₀ = 1.2 and σ₈₅₀ = 0.2 mJy beam⁻¹ with SCUBA-2 in the deep tier of the S2CLS survey, Zavala et al. 2017; σ₈₅₀ = 1.2 mJy beam⁻¹ 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σ, 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 σ_1.1mm = 66.1 μJy beam⁻¹ 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_⊙ yr⁻¹ (SMACS 0723/ALCS) at z = 3–6, obtained by rescaling the 1σ rms with a modified blackbody with temperature T_dust = 40 K, β = 2, k₀ = 0.43 cm² g⁻¹ at λ₀ = 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 ∼3 arcmin² centered on the cluster. AS2UDS and the ALMA archival observations are pointed and covered an area ∼600 × smaller than the parent S2CLS survey in the UDS field (0.96 deg²; Stach et al. 2019). As mentioned in Section 3.4, we retrieve one ∼5σ detection at 850 μm from SCUBA-2 at a 09 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.05}^{+0.09}$); EGS-40032: (3.219, ${3.35}_{-0.11}^{+0.09});$ EGS-31322: (∼3.434, ${3.54}_{-0.10}^{+0.09}$)). 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).

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 (σ_V−J = 0.96 mag). The overlap is less extended when imposing P _{Q ,50%} ≥ 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 05 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.

Zoom In Zoom Out Reset image size

In addition, our selections do not retrieve the dusty candidate QGs at z_phot ∼ 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}_{\mathrm{phot}}={3.58}_{-0.24}^{+0.60}$ and assign it values of ${M}_{\star }={3.0}_{-0.8}^{+1.4}\times {10}^{10}\,{M}_{\odot }$ and P _{Q ,50%} ≪ 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.