RUBIES Spectroscopically Confirms the High Number Density of Quiescent Galaxies from 2 < z< 5

The JWST/NIRSpec Spectroscopy used in this work is taken from the RUBIES program (GO#4233, PIs: A. de Graaff and G. Brammer; A. de Graaff et al. 2025a). RUBIES targets two extragalactic legacy fields: the Extended Groth Strip (EGS) and Ultradeep Survey (UDS), which are parts of the Cosmic Assembly NIR Deep Extragalactic Legacy Survey (N. A. Grogin et al. 2011; A. M. Koekemoer et al. 2011). The parent photometric catalog from which the RUBIES targets were selected was created using public JWST/NIRCam mosaic imaging of EGS from the Cosmic Evolution Early Release Science (CEERS) survey (M. B. Bagley et al. 2023; S. L. Finkelstein et al. 2023, 2025) and in UDS from the Public Release Imaging for Extragalactic Research (PRIMER; C. T. Donnan et al. 2024). RUBIES prioritizes sources that are red in color (F150W − F444W > 2), bright (F444W < 27), or have high photometric redshifts (z_phot > 6.5), which were computed with EAZY (G. B. Brammer et al. 2008). The parent photometric catalog of RUBIES was derived by merging public catalogs that are available on the DAWN JWST Archive (DJA) and was visually inspected to remove invalid entries due to artifacts or pipeline errors. This catalog is used for the completeness correction described in Section 4.3.

All RUBIES spectra used in this work were reduced with msaexp²³ (G. Brammer 2023a), corresponding to version 3 of NIRSpec data (described in detail in A. de Graaff et al. 2025a) released on DJA²⁴ (A. de Graaff et al. 2024; K. E. Heintz et al. 2025). In brief, uncalibrated exposures were taken from the Mikulski Archive for Space Telescopes (MAST)²⁵ and processed through the Detector1Pipeline steps of the standard JWST pipeline, where a mask was inserted for large cosmic-ray snowball events (J. Rigby et al. 2023) calculated with snowblind (J. Davies et al. 2024) and a 1/f correction was applied. After individual slits are identified, the 2D unrectified spectra are flat-fielded and flux calibrated. Generally, RUBIES spectra are reduced with two different background subtraction strategies. However, not all RUBIES spectra have robust spectroscopic redshift quality (grade 3 as defined in A. de Graaff et al. 2025a) in both reduction versions. In one version, the sky background is removed “locally” by taking image differences of the 2D spectra taken at the three spacecraft nod offset positions. Additionally, a second version of RUBIES spectrum reduction was implemented with a global sky subtraction, where a global background is obtained by interpolating spectra from empty slits over the survey footprint and then subtracted from each source spectrum given its spatial location. The global background subtraction version is recommended for bright extended sources. The massive quiescent galaxies that this project aims to search for are bright and expected to be resolved. For any RUBIES spectra included in this analysis (see Section 3.1 for selection criteria), we adopt the global subtraction where available (535 sources), defaulting to local background subtraction for eight sources. Once background subtracted, 1D spectra were optimally extracted (K. Horne 1986) from the rectified 2D spectra (detailed description is available in K. E. Heintz et al. 2025). Finally, using the a priori source position within the shutter and assuming an azimuthally symmetric Gaussian profile, an effective extended-source path-loss correction for light outside the slitlet for each source is derived and applied to the spectrum.

Although this work focuses on the NIRSpec/PRISM spectra as they provide homogeneous coverage of spectral information that constrains the stellar population modeling at 2 < z < 5, we note that RUBIES also collected NIRSpec G395M spectra (2.9–5.3 μm; R ∼ 1000) for these sources. The inclusion of G395M spectra in stellar population synthesis modeling is left for future studies.

The RUBIES photometric targeting catalog was not measured from point-spread function (PSF)-matched NIRCam images and is therefore not optimal for further scientific analysis. For further analysis, we utilize public PSF-matched photometry catalogs of EGS and UDS.²⁶ These catalogs were created with NIRCam F115W, F150W, F200W, F277W, F356W, F410M, and F444W mosaic images of EGS and UDS fields from CEERS and PRIMER (for UDS/PRIMER, the catalog also includes F090W) that are reduced with grizli (G. Brammer 2023b), corresponding to version 7.2 on DJA (F. Valentino et al. 2023). These mosaic images have pixel scales of 0$\mathop{.}\limits^{\unicode{x02033}}$04 pixel⁻¹. All sources were first detected from a long-wavelength sky-subtracted noise-equalized F277W + F356W + F444W stack image, using the Python library for Source Extraction and Photometry (SEP; K. Barbary 2016). Then, PSFs were empirically built in each band, following R. E. Skelton et al. (2014) and K. E. Whitaker et al. (2019), and PSF-matching kernels were produced using Pypher (A. Boucaud et al. 2016) and convolved to degrade all NIRCam images to F444W resolution. Finally, aperture photometry was extracted using circular apertures with SEP (E. Bertin & S. Arnouts 1996; K. Barbary 2016). These aperture fluxes were corrected to total based on each object’s circularized Kron radius and PSF size. A detailed description of photometric catalog construction can be found in J. R. Weaver et al. (2024a). We note that the public versions of the PSF-matched catalogs do not explicitly include rest-frame colors, which we calculate from the eazy-py files using the fsps templates.

The light distributions of all candidate sources exceed the aperture size of JWST’s MSAs. Any potential color gradients can lead to an intrinsic mismatch between the color information obtained from spectroscopy and the standard circular aperture photometry. The spectra mostly target the galaxy cores, while the aperture photometry can include a larger fraction of the full source. To mitigate this potential systematic, we extract a second set of photometry inspired by T. Nanayakkara et al. (2025). This customized photometry approximates the rectangular aperture of the central MSA and is measured from the same PSF-matched 40 mas EGS/CEERS or UDS/PRIMER mosaics. This conservative choice is made because the spectroscopic optimal extraction mostly weights light gathered in the central microshutter. Therefore, we assume that this photometry will be sufficiently close to the region probed by the MSA spectroscopy.

The customized NIRCam photometry further anchors the color information of these systems during modeling when used jointly with the NIRSpec PRISM spectra, which already have outstanding flux calibration. However, the derived luminosity-dependent physical properties, such as stellar mass and star formation rate, would be underestimated because the slitlike aperture cannot capture the total light of these extended objects. Therefore, we take the NIRCam F444W photometry documented in the existing public catalogs (produced with the methodology detailed in J. R. Weaver et al. 2024a) and apply corrections to any spectrophotometric modeling results that depend on the normalization of the SED. We defer the details of such corrections to Section 3.2. For a given source, the public catalogs provide a series of photometry measurements extracted from apertures of various sizes. The catalog also provides a “best-use” value extracted from the recommended aperture size, which varies from source to source depending on the morphology. We default to the photometry from the recommended aperture (i.e., values reported in the SUPER catalog) for the total F444W photometry with which we derive the scaling factor for our models. For sources in this sample, these apertures are typically 0$\mathop{.}\limits^{\unicode{x02033}}$70 or 1$\mathop{.}\limits^{\unicode{x02033}}$00 in diameter.

Although two sets of rest-frame colors (e.g., U – V and V – J) have been used to identify Balmer or 4000 Å breaks in quiescent galaxies, the full SEDs probed by the 0.6–5.3 μm NIRSpec PRISM contain significantly more information. Historically, PCA has been applied to both spectroscopic and photometric data to classify the SEDs of extragalactic sources (e.g., V. Wild et al. 2014). The PCA primarily identifies the main characteristics shared among the dataset, which are known as the eigenspectra. For each galaxy, PCA returns a set of supercolors (SCs), whose total number equals the number of eigenspectra. These SCs are the normalizing coefficients of eigenspectra in the reconstruction of the original spectrum:

$$\begin{eqnarray}&&{F}_{\lambda ,\mathrm{reconstructed}}(\lambda )={F}_{\lambda ,\mathrm{mean}}(\lambda )+\displaystyle \sum _{i=0}^{n-1}{\mathrm{SC}}_{i}\ast {F}_{\lambda ,\mathrm{eigen},i}(\lambda ),\end{eqnarray} \tag{ 1 }$$

where F_{λ,reconstructed}(λ) is the reconstructed spectrum, F_λ,mean(λ) is the mean spectrum, F_λ,eigen,i(λ) are the eigenspectra, and n is the total number of SCs or eigenspectra. As n increases, the reconstructed spectrum converges to the original spectrum. Once a number n is chosen such that the residual between the reconstructed spectrum and the original spectrum is negligible, each set of n SCs can be taken as the low-dimensional representation of each spectrum.

Compared with traditional color selection, PCA enables visualization and classification of observed SEDs without relying on successful fits of spectral synthesis models. At z > 3, early JWST surveys have uncovered novel extragalactic sources that are difficult to constrain by traditional spectral synthesis models, such as “Little Red Dots” (LRDs; e.g., J. E. Greene et al. 2024; J. Matthee et al. 2024). In addition, the classical color selection of quiescent galaxies at this epoch is shown to be incomplete and impure (e.g., S. Belli et al. 2019; J. Antwi-Danso et al. 2023; M. Park et al. 2024; W. M. Baker et al. 2025a). Despite efforts to modify these color criteria, their selection effects are yet poorly understood. PCA provides a data-driven and model-independent approach for selecting and studying the properties of these sources. In this section, we explore the classification of RUBIES SEDs at 2 < z < 5 with PCA. Using PCA classification, we then focus on investigating the preselection of quiescent galaxies that enable the cost-efficient identification of this population.

Although we limit our identification of quiescent galaxies to the RUBIES dataset, for which the targeting strategy is well known, we define the eigenspectra and SC space by leveraging the full set of public NIRSpec/PRISM spectra from the DJA as follows. For RUBIES spectra, we include all spectra that correspond to sources with m_F444W < 24. This ensures a complete selection below this magnitude limit in RUBIES. We further include any PRISM spectra logged with a signal-to-noise ratio (SNR) ≥ 5 on DJA from other JWST spectroscopic programs (including but not limited to JADES, F. D’Eugenio et al. 2025; CEERS, S. L. Finkelstein et al. 2023; UNCOVER, S. H. Price et al. 2025; GTO WIDE, M. V. Maseda et al. 2024).²⁷ All selected spectra have reduction version 3 on DJA and robust spectroscopic redshifts. At low resolution (R ∼ 100 for the NIRSpec PRISM), the selection of quiescent galaxies relies heavily on the information from the Balmer or 4000 Å break as well as the UV slope. Given the wavelength coverage of the NIRSpec/PRISM (0.6–5.3 μm), we adopt a lower redshift cut of z ≥ 2.

We de-redshift and resample these spectra onto a universal rest-frame wavelength grid, using a package implemented in Python, SpectRes (A. C. Carnall 2017). We exclude spectra with significant data discontinuities between 2500 Å < λ_rest-frame < 7400 Å due to reduction quality issues, resulting in a final number of 3206 PRISM spectra. To minimize the contamination of strong emission line flux in our continuum-based selection, we mask out rest-frame 4775–5125 Å ([O III] + Hβ) and 6375–6725 Å ([N II] + Hα + [S II]) in all spectra. We opt to not mask the [O II] λλ3726, 3729 doublet in these spectra. While in some nonquiescent cases there is prominent [O II] emission, the total [O II] equivalent widths are typically smaller than those of [O III] + Hβ or [N II] + Hα + [S II]. Furthermore, due to the coarse wavelength resolution of PRISM around rest-frame 3700 Å, masking this doublet would have reduced spectral sampling of the Balmer break, which is crucial to the selection of quiescent galaxies.

We apply the SparsePCA method implemented by the standard Python package Scikit-learn to our preprocessed spectra. To select the number of eigenspectra n that is the most suitable for our goal of representing and classifying SEDs, we attempt PCA on the same dataset with various n setups. To quantify the quality of the spectroscopic reconstruction, we then define the mean absolute fractional residual of each spectrum as 〈∣(F_λ,original(λ) − F_{λ,reconstructed}(λ))/F_λ,original(λ)∣〉, where the bracket means the average over all wavelength bins. We examine the distribution of the mean fractional residual as a function of n. At n = 4, 99.9% of the sources have a mean fractional residual less than 0.02, suggesting that the reconstructed spectra of these sources account for ∼98% of the flux in the original spectra. At n > 4, each reconstructed spectrum accounts for a higher fraction of flux in the original spectrum, although this improvement becomes marginally better as n grows. At the same time, a higher number of SCs complicates the interpretation of these sources. Therefore, we choose n = 4 for our subsequent analysis.

The choice of SparsePCA is intended to identify eigencomponents corresponding to orthogonal features of the spectral continuum. In total, we identify four eigenspectra that correspond to the UV continuum slope, Balmer or 4000 Å break strength, red optical continuum slope, and IR continuum slope (see Figure 1 for details). This choice allows us to intuitively connect the distribution of the galaxy population in SC space to their physical properties. Thus, we can draw empirical cuts in the SC space to separate galaxies of different spectral types.

Zoom In Zoom Out Reset image size

In Figure 1, we showcase our ability to select objects of different spectral types based on their location in the SC space. The vast majority of the galaxies targeted by early public JWST/NIRSpec programs, and therefore in our PCA, are star-forming galaxies (blue and purple). These sources have positive SC₀, indicating that they possess excessive rest-frame UV flux on top of the mean spectrum. A fraction of them have a notable Balmer or 4000 Å break (blue), as suggested by their slightly negative SC₂. These galaxies have relatively evolved stellar populations and likely had stochastic SFHs (V. Strait et al. 2023; T. J. Looser et al. 2024; A. Covelo-Paz et al. 2025; G. Khullar et al. 2025, in preparation). The LRDs (yellow), whose physical nature is hotly debated, reside in an isolated regime in the SC space. Since they are typically characterized by their unique “V-shaped” continua that turn over at the Balmer limit (e.g., D. J. Setton et al. 2024a; R. E. Hviding et al. 2025), these sources tend to have both positive SC₀ and SC₁. However, the UV colors of faint LRDs are not uniform and consequently, LRDs can shift toward lower SC₀ (C. C. Williams et al. 2024). The quiescent galaxies (green and red) have the strongest Balmer or 4000 Å break strengths and almost no rest-frame UV emission, residing in the lower-left corner of both SC₀–SC₁ and SC₂–SC₃ SC space. The spectral signature of older quiescent galaxies (red) is dominated by stars with longer lifetimes than A-type stars, as they have typically quenched over ∼0.8 Gyr. In contrast, the spectral signature of young quiescent galaxies (or recently quenched and/or poststarburst galaxies; green) is dominated by A-type stars, since they typically quenched within ∼0.8 Gyr. Older quiescent galaxies tend to have flatter optical–NIR continuum slopes and therefore have slightly higher SC₃ and SC₁ than the young quiescent galaxies. They are also characterized by a 4000 Å break instead of the Balmer break characteristic of younger quiescent galaxies, which is reflected by their lower SC₂. However, the exact boundary between dusty star-forming galaxies (maroon) and old quiescent galaxies (red) is difficult to determine in the SC space. As the PCA is essentially an operation that lowers the dimension of the dataset, the subtle difference in the spectral shape required to break the dust–age degeneracy becomes hard to map to the SC location.

Using the visualization of spectral types in Figure 1, we can immediately eliminate sources that are unlikely to be quiescent galaxies without having to perform computationally expensive stellar population synthesis modeling on the entire RUBIES dataset. In order to ensure a complete selection of quiescent galaxies in RUBIES, we adopt generous initial SC cuts, outlined by light-red dashed lines in Figure 1. These cuts effectively remove any sources with spectral features that are mutually exclusive with quiescence. The cut in SC₀ removes sources with significant UV flux, which still host massive stars and therefore recent star formation. The cut in SC₂ removes sources without a prominent Balmer or 4000 Å break, which lack evolved stellar populations. The cuts in SC₁ or SC₃ remove sources with a continuum slope at red optical or NIR wavelengths corresponding to significant dust reddening, which is rare in quiescent galaxies (D. J. Setton et al. 2024; J. C. Siegel et al. 2025).

Using this preselection, we identify an initial sample of 41 unique m_F444W < 24 sources at 2 < z < 5 from RUBIES. Next, we perform spectrophotometric modeling, as follows, to completely disentangle the dusty star-forming and quiescent populations in the initial sample. In Appendix B, we refine these SC cuts to boost selection purity using the best-fitting SED models. The refined SC cuts are shown as red solid lines in Figure 1.

In order to infer the physical properties of PCA-selected galaxies and finalize our selection of quiescent galaxies, we simultaneously fit the PRISM spectrum and NIRCam photometry, using the Bayesian stellar population inference code Prospector (B. Johnson & J. Leja 2017; J. Leja et al. 2017; B. Johnson et al. 2021) with the nested sampling code dynesty (J. S. Speagle 2020). Prospector uses the stellar population synthesis models from the Flexible Stellar Population Synthesis package (C. Conroy et al. 2009; C. Conroy & J. E. Gunn 2010). We adopt the MILES spectral library (P. Sánchez-Blázquez et al. 2006; J. Falcón-Barroso et al. 2011) and MIST isochrones (J. Choi et al. 2016; A. Dotter 2016), assuming a Chabrier IMF (G. Chabrier 2003). Using the ProspectorPolySpecFit instance, we opt to use a polynomial of order 5 to flux calibrate the spectrum to photometry, which we extract from customized apertures described in Section 2.2. We enforce a minimum uncertainty floor of 5% on both the spectrum and photometry.

The setup strategy in our fits is similar to those in A. de Graaff et al. (2025b). We choose a nonparametric SFH that utilizes the continuity prior of Prospector described in J. Leja et al. (2019a). Given the wide range of redshifts in our sample, we adopt different age bins for each object. For the most recent 200 Myr in lookback time, we divide it into four bins of 10, 40, 50, and 100 Myr, then we linearly add four bins of 200 Myr until we reach 1 Gyr in lookback time. The remaining time is evenly divided into N_old bins. We calculate N_old by taking the ceiling of (t_Universe − 1 Gyr)/0.5 Gyr, where t_Universe is the age of the Universe. The number of old bins ranges from one to five for our sample at 2 < z < 5. We fix the redshift to msaexp-derived spectroscopic redshifts. We assume a two-parameter (M. Kriek & C. Conroy 2013) dust law, which is parameterized by the attenuation around old (t > 10 Myr) stars fit in the range τ ∈ [0, 2.5] and a free dust index δ ∈ [−1, 0.4] that describes the deviation from the D. Calzetti et al. (2000) dust law and includes a UV bump that depends on the slope parameterized as in S. Noll et al. (2009). We fix the attenuation around young (t < 10 Myr) stars to be twice that of the older populations. The stellar metallicity is set as a free parameter with a logarithmically sampled uniform prior in the range ${\mathrm{log}}(Z/{Z}_{\odot })\in \,[-2,0.2]$. We mask all wavelengths shorter than rest-frame 1200 Å to avoid contributions from intergalactic medium absorption. We marginalize over the emission lines (specifically, [Ne V] λ3426, [O II] λλ3726, 3729, [O III] λλ4959, 5007, [N II] λλ6548, 6584, [S II] λλ6716, 6731, [S III] λλ9069, 9532, Hα, Hβ, Hγ, and Hδ) by fitting Gaussian profiles. In these Gaussian profiles, the normalization is set as a free parameter, the center is fixed to the line center, and the width is tied to the galaxy intrinsic dispersion and convolved with the instrumental resolution.

The line-spread function curves in the original JWST User Documentation (JDox) are broader than those measured in practice, depending on the exact source morphology (A. de Graaff et al. 2024). Before fitting, all models are convolved with a line-spread function that is a factor of 1.3 narrower than the JDox curves to account for instrumental dispersion, following A. de Graaff et al. (2025b). We additionally include two free velocity dispersions that smooth the stellar continuum and ionized gas emission, which we allow to vary in the range [0, 1000] km s⁻¹ to marginalize over the uncertainty in the line-spread function in addition to the intrinsic dispersion of the galaxy.

As described in Section 2.2, physical properties inferred from SED fitting, such as stellar mass, have to be corrected because the customized photometry does not capture the total light from each galaxy due to the small aperture sizes. We derive an aperture-to-total flux correction from the ratio of F444W flux within the rectangular MSA apertures to the PSF-matched total in the same band. We apply this correction factor to the derived stellar mass and SFH for each galaxy. Once multiplied with the fifth-order polynomial and the scaling factor to F444W total aperture flux, the spectrum in these fits typically increases by a factor of ∼1–2, which is a modest correction since the default slit-loss correction in the DJA pipeline cannot be perfect. Following typical conventions, all masses are reported as surviving stellar mass. We note that this modeling assumes solar-scaled stellar population models. Given that these massive galaxies are probably α-enhanced (e.g., A. G. Beverage et al. 2025), this will slightly overestimate the masses and ages (M. Park et al. 2024).

In Figure 2, we show two examples of our spectrophotometric fits: a young quiescent (poststarburst) galaxy at z ∼ 4 (top) and an old quiescent galaxy at z ∼ 2.7 (bottom). In the upper left panels, the observed spectrum and its uncertainties are shown in red, and the observed photometry with error is shown in orange. The best-fitting model spectrum and photometry are shown in black. Both the observed and model spectra have been scaled by the polynomial calibration vector to match the photometry. The insets show the NIRCam/F444W cutout (2$\mathop{.}\limits^{\unicode{x02033}}$4 in width) of the target galaxy. The red rectangle traces the central MSA (for slitlike aperture photometry), and the red circle traces the circular aperture from which the catalog photometry was extracted. The bottom left panels show the residuals for the spectrum and photometry of the fits (χ values), which we define as the difference between observed flux and best-fitting model flux, normalized by the uncertainties in the observed flux. The right panels show the median and 68% confidence interval SFHs for these galaxies as solid black curves and gray regions, respectively. We show the median and 68% confidence interval of t₉₀ (lookback time at which 90% of the current stellar mass was formed) in the same panels as dashed red lines and light-red regions, respectively. Overall, these fits are excellent and enable the characterization of SFHs and the physical properties of these galaxies. The SED fits and SFHs of the remaining galaxies in this sample are shown in Appendix A.

Zoom In Zoom Out Reset image size

Based on our initial PCA selection, we find that that the majority (∼90%) of the bright (F444W < 24 ∪  SNR > 5) RUBIES sources at 2 < z < 5 possess empirical spectral features that unambiguously eliminate their possibility of being quiescent systems (see Section 3.1 for a qualitative discussion). Our approach to identifying quiescent galaxies is iterative. First, we apply the initial conservative SC cuts (light-red dashed lines in Figure 1) to focus on potentially quiescent sources and minimize the computational cost. We model the preselected galaxies with Prospector and infer their physical properties. We examine the relation between the SCs and specific SFR (sSFR) among these galaxies and identify a smaller region in SC space that robustly identifies quiescent galaxy SEDs (red solid lines in Figure 1). This selection largely avoids contamination from other red sources, including dusty star-forming galaxies and those with decreased, but significant, ongoing star formation.²⁸ Next, we describe the identification and classification of quiescent galaxies, using the physical properties and SFHs of galaxies selected by the refined SC cuts.

In the left panel of Figure 3, we show the best-fitting stellar mass versus sSFR from Prospector for all galaxies in the refined SC-selected sample. To test whether binary classification is sensitive to measurement uncertainties, we define two quiescent criteria as follows:

1. 1.  
   Criterion 1 (red symbols). sSFR_50th < 10⁻¹⁰ yr⁻¹ and
2. 2.  
   Criterion 2 (pink symbols). sSFR_16th < 10⁻¹⁰ yr⁻¹.

We adopt Criterion 1 as our primary criterion of quiescence, which yields 17 quiescent galaxies in total as our fiducial sample. Criterion 2 includes three more galaxies that could have been selected by Criterion 1 if their sSFRs are perturbed by 1σ. In the following discussion, we also refer to these three galaxies as “marginally quiescent.” We verify that an evolving sSFR threshold (e.g., sSFR_50th < 0.2/t_Universe(z) Gyr⁻¹ as in A. C. Carnall et al. 2023a; W. M. Baker et al. 2025a) identifies almost the same quiescent galaxies as our fiducial sample, including one additional galaxy at 4 < z < 5 (ID: RUBIES-UDS-140707).

Zoom In Zoom Out Reset image size

Using the SFHs inferred by Prospector, we calculate the formation timescale (t₉₀ or the lookback time by which a galaxy formed 90% of its current stellar mass) of each quiescent galaxy. In the following analysis, we use t₉₀ as a proxy for time since quenching to differentiate between the old and young galaxy populations. In the middle panel of Figure 3, we show t₉₀ versus stellar mass for all quiescent galaxies. Most of the quiescent galaxies in our sample are young (pink or red) with best-fit t₉₀ < 0.8 Gyr. In the right panel of Figure 3, we show the distribution in stellar mass versus redshift. All older (t₉₀ > 0.8 Gyr) sources (black symbols) are found at the lowest redshifts (2 < z < 3). The SEDs of these galaxies prominently display 4000 Å breaks rather than Balmer breaks (see Figures 2 and 8). The earliest emergence of the older quiescent population at z ∼ 3 in this sample is roughly contemporaneous with the record holder of the “maximally old” quiescent galaxy at z = 3.2 reported in K. Glazebrook et al. (2024). At face value, they are likely the descendants of the first generation of massive quiescent galaxies at z ≥ 5–7 (A. de Graaff et al. 2025b; A. Weibel et al. 2025). However, we note that the stellar masses of these examples are slightly lower than the current higher-redshift counterparts in A. Weibel et al. (2025) and A. de Graaff et al. (2025b). This discrepancy could simply reflect cosmic variance or a systematic bias in mass inference between galaxies in different evolution stages. The inferred physical properties relevant to the selection and classification of the complete quiescent sample are tabulated in Table 2.

This final sample includes the most massive quiescent galaxies (median M_* ∼ 10^10.7 M_⊙) at 2 < z < 5 in RUBIES, in part because we enforce the parent sample to include only the bright galaxies with m_F444W < 24. However, this magnitude limit does not translate to a uniform mass limit for the entire final sample. Although these galaxies span a wide range of redshifts (2 < z < 5), galaxies are younger and therefore have higher mass-to-light ratios at higher redshift, somewhat mitigating the increasing luminosity distance.

To derive the redshift-dependent effective mass–completeness limit, we obtain rescaled analogs of the observed galaxies, using a grid of redshifts and stellar masses and best-fitting Prospector models. For a group of models within a range of redshift, we compute their F444W magnitude corresponding to a given stellar mass and examine where they would pass the F444W < 24 threshold. The details of this analysis are presented in Appendix B. It suggests that our magnitude-limited selection is complete above 10^10.3 M_⊙ at 2 < z < 3, above 10^10.5 M_⊙ at 3 < z < 4, and above 10^10.6 M_⊙ at 4 < z < 5.

We also investigate how dust attenuation (A_V) could affect the locations of the true quiescent population in this sample in SC space. Using the best-fit Prospector models for all sSFR_50th < 10⁻¹⁰ yr⁻¹ galaxies, we re-generate model SEDs for a grid of A_V values. The details of this analysis are presented in Appendix B. In brief, the true quiescent galaxies generally move along a diagonal track with a negative slope in SC₀–SC₁ and SC₂–SC₃ as their dust attenuation increases. In addition to the initial SC thresholds, the diagonal SC cuts in SC₀–SC₁ and SC₂–SC₃ (shown as red solid lines in Figure 1) are roughly parallel to the A_V evolution track of the true quiescent population. We find that our refined SC cuts should be complete for quiescent galaxies with A_V < 0.7.

While spectroscopy remains the most robust tool to identify quiescent galaxies, photometric samples will always provide a comprehensive and cost-efficient probe of the cosmic volume, although at the expense of purity (e.g., A. C. Carnall et al. 2023a; F. Valentino et al. 2023; S. Alberts et al. 2024; A. S. Long et al. 2024). We now estimate the purity and completeness of various rest-frame color selection criteria if applied to this spectroscopic sample of massive quiescent galaxies at 2 < z < 5.

In Figure 4, we show the rest-frame UVJ color distribution of all galaxies in EGS and UDS with F444W < 24 at 2 < z < 5 in the left panel and their rest-frame u_sg_si_s colors (J. Antwi-Danso et al. 2023) in the middle panel. The galaxies in RUBIES that are selected by our refined SC cuts are shown as hexagon symbols, whose color coding follows previous convention. Other galaxies in RUBIES that are excluded by the refined SC cuts are shown as black dots. All the remaining field galaxies not surveyed by RUBIES are binned into 2D histograms, where the darker color indicates the higher density of objects.

Zoom In Zoom Out Reset image size

The conventional rest-frame UVJ selection criteria with a horizontal U – V cut (e.g., R. J. Williams et al. 2009; K. E. Whitaker et al. 2011; C. Schreiber et al. 2015) misses some of the youngest quiescent galaxies in our sample, as expected (e.g., W. M. Baker et al. 2025a). Extending to bluer U − V colors (diagonal dashed lines; S. Belli et al. 2019; W. M. Baker et al. 2025a) captures many excluded galaxies, for this quiescent sample (sSFR_50th < 10⁻¹⁰ yr⁻¹) ∼30% have U − V < 1.23. We find contamination fractions of ∼35% in both the classical and extended UVJ criteria, which is high but largely consistent with the ∼10%–30% reported in J. Leja et al. (2019b), J. Antwi-Danso et al. (2023), and T. Nanayakkara et al. (2025). We return to discuss the impact of false positivity on number densities from UVJ color selections in Section 6. Alternatively, the u_sg_si_s criteria select almost all (∼95%) quiescent galaxies in this sample at the cost of much higher contamination rates (∼60%).

Even though the photometric redshifts of these bright quiescent galaxies agree well with the spectroscopic redshifts (see right panel of Figure 4), the substantial uncertainties in photometric redshifts can contribute significant systematic uncertainties to the rest-frame color estimations and scatter in color–color space. The typical uncertainties in these color estimations are shown in the left and middle panels of Figure 4. We emphasize that this scatter is not included in our purity and completeness fractions quoted above. Both the impostors and true quiescent galaxies located near these selection boundaries can shift into or out of these color selection criteria if perturbed within the confidence interval of their rest-frame colors. These uncertainties and the nuance in the definition of true quiescence (Criterion 1 versus Criterion 2) complicate the determination of systematic impurity or incompleteness in these rest-frame color selections.

Another caveat is that the completeness or purity fraction reported in this section only considers the EGS and UDS galaxies surveyed by RUBIES and is not corrected for the selection bias introduced by the RUBIES selection function. Without this correction, the purity fraction of the extended UVJ criteria is likely prone to underestimation. These criteria include younger quiescent galaxies but also bluer impostors (K. E. Whitaker et al. 2012; L. R. Spitler et al. 2014; C. M. S. Straatman et al. 2015). Since RUBIES preferentially targeted redder sources, which we discuss in detail in the following section, RUBIES could miss more blue impostors and underestimate the total number of sources present in the extended color region.

Spectroscopy reliably prevents the contamination of dusty star-forming galaxies in the selection of quiescent galaxies, which is unavoidable with only photometry. However, these spectroscopic surveys cannot comprehensively target every eligible source in the field. The number densities derived from this kind of spectroscopic selection must be statistically corrected for sample incompleteness. For the RUBIES massive quiescent galaxy sample, the overall sample completeness is a net product of three selection steps.

First, we consider the fraction of field sources that are eligible to be included as survey candidates (S_el). This is a function of the apparent photometric properties of the source ${\boldsymbol{m}}=({m}_{{\rm{F}}444{\rm{W}}},{m}_{{\rm{F}}150{\rm{W}}}-{m}_{{\rm{F}}444{\rm{W}}})$, which depends on their intrinsic physical properties such as stellar mass (M_*), sSFR, redshift (z), and dust reddening (A_V). Second, we define a variable, S_slit, to be the fraction of potential sources for which suitable spectra were obtained. This fraction includes whether a source is assigned an MSA slit, whether a spectrum was successfully obtained and reduced, and whether the spectroscopic wavelength coverage, redshift quality grade, and SNR meet the PCA requirements described in Section 3.1. This is similarly a function of the photometric properties of the source ${\boldsymbol{m}}$. Third, we consider the fraction of sources included in our SC cuts among all potential quiescent sources for which a successful spectrum is obtained (S_PCA). This term is a function of the source spectrum f_λ and accounts for missing potential quiescent galaxies with unusual SCs. In summary, we express these terms as:

1. (a)  
   ${S}_{\mathrm{el}}({\boldsymbol{m}}({M}_{\ast },\mathrm{sSFR},z,{A}_{{V}}))$,
2. (b)  
   ${S}_{{\rm{slit}}}({\boldsymbol{m}}({M}_{\ast },{\rm{sSFR}},z,{A}_{V}))$, and
3. (c)  
   S_PCA(f_λ(M_*, sSFR, z, A_V)).

Therefore, the total sample completeness is given by

$$\begin{eqnarray}&&{S}_{\mathrm{tot}}({M}_{* },\mathrm{sSFR},z,{A}_{V})={S}_{\mathrm{el}}\cdot {S}_{\mathrm{slit}}\cdot {S}_{\mathrm{PCA}}.\end{eqnarray} \tag{ 2 }$$

We simplify some of the terms in the selection function as follows. Since RUBIES considers any field source with F444W < 28.5 to be eligible and we have enforced all sources in this sample to have m_F444W < 24, we assume S_el = 1 for all galaxies in this sample. Furthermore, S_PCA has no dependency on stellar mass or redshift, since our PCA is conducted on normalized, rest-frame spectra. As we have shown, the SC selection of quiescent galaxies has an effective limit in dust attenuation of A_V < 0.7, where we assume S_PCA ∼ 1 at A_V < 0.7 and sSFR < 10⁻¹⁰ yr⁻¹. With these simplifications, we show that S_tot only depends on S_slit for the quiescent galaxies in our selection.

RUBIES selects spectroscopic targets from eligible sources in CEERS/EGS and PRIMER/UDS according to their photometric redshifts, F444W magnitudes, and F150W − F444W colors (A. de Graaff et al. 2025a). To estimate S_slit, we divide the finalized quiescent galaxies, RUBIES sources with successful spectra, and eligible parent sources into three redshift bins: [2, 3], [3, 4], and [4, 5]. In Figure 5, we show the color–magnitude distribution of finalized quiescent galaxies (using the same color coding as in previous figures), RUBIES sources with successful spectra (black dots), and the parent photometric catalog sources (gray dots) in each redshift bin. We select the approximate color–magnitude space occupied by the quiescent galaxies in this sample in each redshift bin (thick gray dashed line) and divide these color–magnitude spaces into smaller boxes that are 1 mag in width (thin gray dashed line). We mask out the faint regime (m_F444W > 24), which is excluded by the quiescent galaxy selection in this sample. In each color–magnitude box (m), we count the total number of quiescent galaxies in this sample N_{QG,surveyed,m} (hexagons), the number of RUBIES targets with a successful spectrum (N_surveyed,m; black dots), and the number of eligible parent sources (N_total,m; gray dots). The fraction of eligible sources which yield successful spectra varies smoothly in the color–magnitude space (A. de Graaff et al. 2025a). This allows us to approximate this fraction as a constant per color–magnitude box: ${S}_{{\rm{slit}}}({\boldsymbol{m}})\sim {N}_{{\rm{surveyed}},{\boldsymbol{m}}}/{N}_{{\rm{total}},m}$, when ${\boldsymbol{m}}\in {\boldsymbol{m}}$.

Zoom In Zoom Out Reset image size

A few sources selected by RUBIES to be assigned an MSA slit fail to produce a spectrum that meets the PCA prerequisite in wavelength coverage, redshift quality grade, or SNR. However, we find that this failure rate does not have a significant dependence on color, magnitude, or redshift within the relevant parameter space. Generally, the spectroscopic redshifts recovered by RUBIES in the entire survey agree well with photometric redshifts (the normalized median absolute deviation between the two is small: σ(Δz/(1 + z)) = 0.033; A. de Graaff et al. 2025a). For the bright quiescent galaxies in this sample, the match between photometric and spectroscopic redshifts is excellent, as shown in the right panel of Figure 4. We assume that any remaining unsurveyed quiescent galaxies should also have reliable photometric redshifts. Therefore, the approximation described above is likely not affected by any catastrophic errors in photometric redshifts or biased by any potential correlation between spectrum acquisition failure and source color, magnitude, or redshift.

For a population of RUBIES massive quiescent galaxies in a narrow redshift bin, which have ${\rm{sSFR}}\lt {{\rm{sSFR}}}_{{\rm{\max }}}$, ${M}_{* }\gt {M}_{* ,\min }$, and ${A}_{{V}}\lt {A}_{V,\max }$, the number density is given by their expected number divided by the effective survey volume:

$$\begin{eqnarray}&&\begin{array}{l}n({M}_{* }\gt {M}_{* ,\min },\mathrm{sSFR}\lt {\mathrm{sSFR}}_{\max },{A}_{{V}}\lt {A}_{V,\max })\\ =\,\displaystyle \frac{1}{{V}_{\mathrm{eff}}}{\int }_{{M}_{* ,\min }}^{\infty }{\int }_{-\infty }^{{\mathrm{sSFR}}_{\max }}{\int }_{0}^{{A}_{V,\max }}\\ \displaystyle \frac{{\rm{\Phi }}({M}_{* },\mathrm{sSFR},{A}_{{V}})}{{S}_{\mathrm{tot}}({M}_{* },\mathrm{sSFR},{A}_{{V}})}d{M}_{* }\,d\mathrm{sSFR}\,d{A}_{{V}},\end{array}\end{eqnarray} \tag{ 3 }$$

where Φ(M_*, sSFR, A_V) is the number of RUBIES galaxies per stellar mass per sSFR per attenuation per comoving cosmic volume. We do not intend to solve for this integral in the intrinsic parameter space, since the true distribution of Φ(M_*, sSFR, A_V) is challenging to constrain. We instead assume that the above integral in the intrinsic parameter space (within the effective selection limits) is equivalent to the following summation in the observed parameter space:

$$\begin{eqnarray}&&{n}_{\mathrm{Massive},Q}=\displaystyle \frac{1}{{V}_{\mathrm{eff}}}\sum ^{{\boldsymbol{m}}}\displaystyle \frac{{N}_{\mathrm{QG},\mathrm{surveyed},{\boldsymbol{m}}}\cdot {N}_{\mathrm{total},{\boldsymbol{m}}}}{{N}_{\mathrm{surveyed},{\boldsymbol{m}}}},\end{eqnarray} \tag{ 4 }$$

where n_Massive,Q is the number density of massive quiescent galaxies. And V_eff is given by

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}=\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{field}}}{3}\cdot [{d}_{\mathrm{com}}{({z}_{\max })}^{3}-{d}_{\mathrm{com}}{({z}_{\min })}^{3}],\end{eqnarray} \tag{ 5 }$$

where Ω_field is the angular size of the effective survey area and d_com(z) is the comoving distance to redshift z.

This approximation assumes that any source with similar intrinsic properties (M_*, sSFR, and A_V) should display similar observed source properties (m_F444W and m_F150W − m_F444W), in the special case of massive (M_* > 10^10.3 M_⊙) quiescent (sSFR < 10⁻¹⁰ yr⁻¹) galaxies without significant dust reddening (A_V < 0.7). We essentially estimate how the selection of quiescent galaxies in RUBIES depends on their intrinsic source properties, given that we know how the selection of RUBIES sources depends on their observed source properties.

We show the resulting comoving number densities in the left panel of Figure 6 in three flavors. The RUBIES survey footprints overlap half of the sky covered by JWST NIRCam imaging in EGS and UDS (A. de Graaff et al. 2025a). In part to maximize survey efficiency, the MSA mask design of RUBIES targets overdensities within these fields and could result in biases in our number density estimations. To probe how potential cosmic overdensities under the RUBIES footprint can impact our conclusions, we take the fiducial quiescent sample (sSFR_50th < 10⁻¹⁰ yr⁻¹) and calculate the number density in each redshift range using two different approaches. In the first approach, we only count parent catalog sources that fall under the survey footprint for completeness estimation. In this case, we divide the resulting expected number of quiescent galaxies by the comoving volume corresponding to only the total sky area under the footprints (${{\rm{\Omega }}}_{{\rm{field}}}\sim 150\,{{\rm{arcmin}}}^{2}$). In the second approach, we count sources from the entire UDS/EGS footprints as N_total for completeness estimation. In this case, we use the comoving volume corresponding to the entire sky area of the NIRCam imaging (${{\rm{\Omega }}}_{{\rm{field}}}\sim 304\,{{\rm{arcmin}}}^{2}$). As shown in the left panel of Figure 6, the targeting bias due to overdensities at 4 < z < 5 results in a footprint-only number density (red) ∼0.2 dex higher than that of the full field (orange). This offset effectively reflects the cosmic variance within EGS and UDS.

Zoom In Zoom Out Reset image size

To test the impact of different quiescent criteria in Section 4.1, we calculate a third flavor with the extended quiescent sample (sSFR_16th < 10⁻¹⁰ yr⁻¹), adopting the “footprint-only” approach. An extended criterion for quiescent galaxy selection (Criterion 2) slightly increases the resulting number density (pink; Figure 6, left panel) at all redshifts by ∼0.1 dex. We tabulate the various flavors of quiescent galaxy number densities in Table 1. Overall, these variations are insignificant compared to uncertainties due to cosmic variance, which we will describe in the following. We choose to report the number densities calculated with the “Criterion 1” sample and the “footprint-only” approach as our fiducial result.

Table 1. RUBIES Massive Quiescent Galaxy Number Densities

Mass or Magnitude Limit	mF444W < 24	mF444W < 24	M* > 1010.5 M⊙
Quiescence Criterion	sSFR50th < 10−10 yr−1	sSFR16th < 10−10 yr−1	sSFR50th < 10−10 yr−1
	(Mpc–3)	(Mpc–3)	(Mpc–3)
2 < z < 3	$1.1{0}_{-0.24}^{+0.56}\cdot 1{0}^{-4}$	$1.2{0}_{-0.26}^{+0.61}\cdot 1{0}^{-4}$	$0.7{7}_{-0.19}^{+0.48}\cdot 1{0}^{-4}$
3 < z < 4	$3.5{3}_{-0.57}^{+0.97}\cdot 1{0}^{-5}$	$4.0{5}_{-0.68}^{+1.17}\cdot 1{0}^{-5}$	$2.7{4}_{-0.48}^{+0.86}\cdot 1{0}^{-5}$
4 < z < 5	$1.2{6}_{-0.38}^{+1.01}\cdot 1{0}^{-5}$	$1.5{7}_{-0.47}^{+1.21}\cdot 1{0}^{-5}$	$1.2{6}_{-0.38}^{+1.01}\cdot 1{0}^{-5}$

Download table as:  ASCIITypeset image

To estimate the measurement uncertainties on the derived number densities, we calculate the binomial proportion confidence interval for the spectroscopic completeness in each color–magnitude box (N_surveyed,m/N_total,m), following the method in E. B. Wilson (1927; the Wilson score interval). We propagate the confidence interval in completeness into number densities, using Equation (4). We plot the resulting uncertainties due to counting and completeness corrections as colored error bars in the left panel of Figure 6.

Additionally, we estimate the systematic uncertainty due to cosmic variance using a method similar to that described in L. Taylor et al. (2022). Briefly, we use ∼1300–2300 mock light cones of similar area to the RUBIES survey derived from the UniverseMachine model (P. Behroozi et al. 2019) run on the MDPL2 dark-matter-only N-body simulations (K. Riebe et al. 2013). Mock galaxies are rank ordered by stellar mass, and comparison samples are abundance matched (to observed number densities in each redshift range) to estimate the 1σ scatter in the field-to-field variation. Due to the small survey volume of RUBIES, the estimated cosmic variance is large (∼60%) and dominates the uncertainty budget for this sample.

In the middle panel of Figure 6, we compare our fiducial quiescent galaxy number densities to other JWST observations. Overall, our measurements are in good agreement with other JWST results. Notably, W. M. Baker et al. (2025a), which is the only other spectroscopic sample shown, is the most consistent with our results at 2 < z < 4. That work combined multiple NIRSpec observation programs (D. J. Eisenstein et al. 2025, 2023) targeting the GOODS fields (M. Giavalisco et al. 2004) and utilizes a similar approach: preselecting quiescent candidates with rest-frame UVJ colors (with an extended color cut) and finalizing with results from spectrophotometric modeling. They report five quiescent galaxies at 4 < z < 4.5 but none above z > 4.5. Taken together with our results, this suggests that the number density of quiescent galaxies likely drops steeply at 4 < z < 5.

The other literature samples shown in the middle panel of Figure 6 are derived from photometric samples. A. C. Carnall et al. (2023a) selected quiescent galaxies based on sSFRs inferred from photometry-only SED fits, using the NIRCam imaging data of CEERS/EGS (which partially overlaps with the RUBIES footprint); we show the “robust, bright” subsample.²⁹ That study adopted a redshift-dependent sSFR cut that roughly approximates a rest-frame UVJ selection. A. S. Long et al. (2024) selected quiescent galaxies from the same dataset (CEERS/EGS) based on observed-frame colors, using color cuts informed by model templates. F. Valentino et al. (2023) derived number densities by combining the quiescent galaxies selected with rest-frame colors from all available fields at the time. We choose the version of their results computed with the massive subsample, which was selected with the padded UVJ selection criteria.

Most of the number densities derived from photometric samples appear to be (0.3–0.5 dex) higher than ours at 3 < z < 5. These studies overlap with RUBIES partially in CEERS/EGS or PRIMER/UDS. While the nonoverlapping sources between RUBIES and these studies can contribute to this discrepancy, it is also likely that the discrepancy between the A. C. Carnall et al. (2023a) sample and RUBIES is driven by the contamination of dusty star-forming objects (∼10%–35% in rest-frame UVJ selection, as reported in J. Leja et al. 2019b; J. Antwi-Danso et al. 2023; T. Nanayakkara et al. 2025; and Section 4.2). We note that A. S. Long et al. (2024) use empirical colors to select quiescent galaxies for which the systematic contamination or incompleteness is challenging to estimate. F. Valentino et al. (2023) report a slightly lower number density at 3 < z < 4, which can be explained by its mass limit being 0.1 dex higher than the effective stellar mass limit in our sample at this redshift. In addition, the known overdensities at 3 < z < 4 in CEERS/EGS (“Cosmic Vine;” S. Jin et al. 2024) can drive the number densities in this work to be higher, since the F. Valentino et al. (2023) sample is further constrained by various other lines of sight. However, the contribution to the final statistic from intrinsic overdensities in CEERS is likely small, since similar results are obtained in W. M. Baker et al. (2025a).

In the right panel of Figure 6, we separate the Criterion 1 galaxies in this sample by t₉₀ into young and old quiescent galaxies and calculate their number densities respectively using parent catalog sources in the RUBIES footprint. This is a key benefit of including spectroscopic information in our analysis. We further compare these to (pre-JWST) literature measurements at lower redshifts that rely on spectroscopy (B. Forrest et al. 2018; S. Belli et al. 2019) or photometry (M. Clausen et al. 2024) to estimate galaxy ages. We find that the number density of older quiescent galaxies from pre-JWST samples is systematically lower than the median estimation from RUBIES at 2 < z < 3, although these results are marginally within our measurement uncertainties. However, for young quiescent galaxies, there is a large discrepancy (almost 1 dex) in number density at z > 3 between RUBIES and the only non-JWST result at these redshifts (B. Forrest et al. 2018), potentially due to differences in effective mass limits. We find a relatively consistent number density (within 0.3 dex) with non-JWST literature results with similar minimum mass at 2 < z < 3 (B. Forrest et al. 2018; M. Clausen et al. 2024). The young quiescent galaxy number density at 2 < z < 3 reported by S. Belli et al. (2019) is much lower (∼0.8 dex), but for a much more massive sample (M_* > 10^10.8 M_⊙ versus M_* > 10^10.3 M_⊙ for RUBIES).

At z ∼ 3, these ground-based samples likely miss both the old and young quiescent galaxies compared to any JWST sample such as RUBIES, due to several factors. The selection of pre-JWST-era samples at these epochs, such as the one in B. Forrest et al. (2018; following M. Kriek et al. 2011), primarily relied on magnitude-limited K-band surveys. At low mass or high redshifts, quiescent galaxies will fall below K-band selection limits. This is especially relevant for the oldest galaxies, which have the highest M/L in the K band. This incompleteness at z > 3 is further exacerbated by the lack of rest-frame NIR coverage, which requires extrapolation to compute rest-frame V − J colors for common UVJ color selection (J. Antwi-Danso et al. 2023). Meanwhile, JWST samples reach deeper magnitude limits, cover longer wavelengths, and have yielded relatively consistent number densities of quiescent galaxies at z > 3.

Constraining the number density of the first old quiescent galaxies at cosmic noon provides an opportunity to validate the number density of young quiescent galaxies at earlier epochs. After approximately ∼0.5–1 Gyr, young quiescent galaxies at a given epoch lose the A-type stars that dominate their SEDs and thus become part of the old quiescent population in the next epoch (K. E. Whitaker et al. 2012). Therefore, the old quiescent population should be a cumulative ensemble of all quiescent galaxies quenched ∼0.5–1 Gyr ago, if rejuvenation or major mergers are infrequent among these galaxies. Since the number density of young quiescent galaxies appears to drop steeply at 4 < z < 5 (as previously discussed) and the Universe is only ∼1 Gyr old at z ∼ 5, we expect the majority of the old quiescent population uncovered at 2 < z < 3 to be descendants of the young population at z ∼ 4. Our measurements are consistent with this picture: the number density of old quiescent galaxies at 2 < z < 3 is indeed similar to or slightly higher than those of the young galaxies ∼1 Gyr beforehand. However, we note that the number density of old quiescent galaxies is fairly uncertain (∼0.5 dex) even with a broad redshift bin (equivalent to ∼1 Gyr). In the future, a larger spectroscopic census of the old quiescent population at 2 < z < 3 will help further test this interpretation.

In Figure 7, we compare the observed number density of quiescent galaxies derived from this sample to various simulation predictions discussed in C. d. P. Lagos et al. (2025). Three are semianalytical models: Shark (C. d. P. Lagos et al. 2018, 2024), GAEA (G. De Lucia et al. 2014, 2024; M. Hirschmann et al. 2016), and Galform (C. G. Lacey et al. 2016). The other three are cosmological hydrodynamical simulations: Eagle (R. A. Crain et al. 2015; J. Schaye et al. 2015; S. McAlpine et al. 2016), IllustrisTNG (A. Pillepich et al. 2018; V. Springel et al. 2018), and Simba (R. Davé et al. 2019). All three semianalytical models correspond to comoving volumes of ∼700 cMpc³, and the three cosmological hydrodynamical simulations have comoving box sizes of ∼100–150 cMpc³. All simulation predictions shown here assume a quiescence criterion of sSFR < 10⁻¹⁰ yr⁻¹ and a stellar mass limit of ${\rm{log}}({M}_{* }/{M}_{\odot })\gt 10.5$, regardless of dust attenuation. The simulation number densities shown in the left panel select galaxies by their exact values in stellar mass and SFR, while those in the right panel additionally consider the random errors in these properties, as described in C. d. P. Lagos et al. (2025) and G. De Lucia et al. (2024). For the curves shown in the right panel, number densities are recalculated after convolving the galaxy mass function or SFR distribution in these simulations with a Gaussian to mimic the scattering of the galaxy population in mass or SFR due to errors. We assume errors in stellar mass and SFR to be independent and these Gaussians are centered at zero with widths of 0.3 dex (except for GAEA, where the widths are 0.35 dex). These choices of Gaussian widths for error convolution are motivated by the typical uncertainties in stellar mass and SFR inferred from multiwavelength observations (A. S. G. Robotham et al. 2020; S. Bellstedt & A. S. G. Robotham 2025). Notably, incorporating scatter to emulate the effects of measurement uncertainties systematically increases the number densities in Eagle and Shark. This is likely because >10^10.5 M_⊙ falls in the exponential decline of the quiescent galaxy stellar mass functions in these simulations (C. d. P. Lagos et al. 2025), introducing a net upward Eddington bias and inflating number densities. For the observed quiescent galaxy number density shown in these panels, we take the fiducial sample (Criterion 1) and remove galaxies with ${\mathrm{log}}({M}_{* }/{M}_{\odot })\lt 10.5$ in the 2 < z < 3 bin, in order to be consistent with the sSFR or mass limit in these simulations. W. M. Baker et al. (2025b) have shown that the observed quiescent galaxy stellar mass function is flat at a ${\mathrm{log}}({M}_{* }/{M}_{\odot })\sim 10.5$ cutoff at 2 < z < 4. Therefore, the number density derived from RUBIES is potentially not sensitive to the Eddington bias discussed above at z < 4.

Zoom In Zoom Out Reset image size

At face value, only Simba, IllustrisTNG, and GAEA agree well with the observed number density of quiescent galaxies at 2 < z < 3. After emulating measurement uncertainties, all six simulations are largely consistent with observational data at cosmic noon. However, at earlier times, all simulations underpredict the observed number densities, by ∼0.4 dex at 3 < z < 4 and ∼1 dex at 4 < z < 5. This discrepancy persists regardless of our empirical definition of quiescence; for example, using the evolving quiescence criterion of sSFR < 0.2/t_Universe(z) Gyr⁻¹ produces similar results. Given the median theoretical number density, we would expect to find ∼1 quiescent galaxy at ${\mathrm{log}}({M}_{* }/{M}_{\odot })\gt 10.5$ at 4 < z < 5 in the total sky area covered by EGS and UDS ($304\,{{\rm{arcmin}}}^{2}$). Yet we find five such galaxies in the ∼50% covered by the RUBIES survey ($150\,{{\rm{arcmin}}}^{2}$). We note that cosmic variance plagues simulations and observations alike. These cosmological simulations have small box sizes (∼100 cMpc³), for which a number density of ∼10⁻⁶ cMpc⁻³ corresponds to one object in the entire simulation box. At 4 < z < 5, their number densities are sensitive to random counting errors. However, we emphasize that this comparison is still meaningful since the number density of quiescent galaxies in the RUBIES sample is an order of magnitude higher (∼10⁻⁵ cMpc⁻³). Given our consistency with previous studies, we conclude that the dramatic discrepancy between the observed and predicted quiescent galaxy populations before z > 3 is unlikely to be attributed to contamination within photometric samples or targeting biases in spectroscopic studies.

Although we expect the effect to be small, our measured number densities could be slightly underestimated at the highest redshifts due to the effective magnitude completeness limits within the RUBIES survey. The magnitude limit (F444W < 24) of our fiducial parent sample corresponds to stellar mass limits of 10^10.3 M_⊙ (2 < z < 3), 10^10.5 M_⊙ (3 < z < 4), and 10^10.6 M_⊙ (4 < z < 5), as discussed in Appendix B. Therefore, RUBIES could have failed to target galaxies $1{0}^{10.5}\lt {\mathrm{log}}{M}_{* }\lt 1{0}^{10.6}\,{M}_{\odot }$ with high intrinsic mass-to-light ratios or near z ∼ 5. We expect this effect to be insignificant given our simulations (Figure 10). We note that it is further possible, though unlikely, that this sample is missing a significant population of heavily dust-attenuated A_V > 0.7 quiescent galaxies (see details in Appendix B). Although most quiescent galaxies at z < 2.5 have A_V < 0.75 (e.g., K. A. Suess et al. 2019; J. C. Siegel et al. 2025), rare counterexamples with significant dust reddening in the core (A_V > 0.75, D. J. Setton et al. 2024; J. C. Siegel et al. 2025) or the outskirts (Z. Ji et al. 2024) exist. However, these sources of sample incompleteness would only exaggerate the tension between simulations and observations.

Among these simulations and models, the wide range of predicted quiescent number densities at a given redshift is mainly due to the different implementations of AGN feedback (see C. d. P. Lagos et al. 2025, for detailed discussion). Further modifications to the AGN feedback implementation could resolve the current tension in quiescent galaxy number densities. Although hard to pinpoint observationally, many lines of evidence point toward the simultaneity of quenching and AGN activity. For example, the AGN incidence rate of massive quiescent galaxies at similar epochs (∼50% from a multiwavelength search by W. M. Baker et al. 2025a; ∼20% from analysis of optical emission lines in M. Martínez-Marín et al. 2024; K. Ito et al. 2025a) is much higher than that of the youngest quiescent galaxies at low redshifts (J. E. Greene et al. 2020). A number of the quiescent galaxies in this sample exhibit strong nebular emission lines, which hints at the incidence of nuclear activities, although we defer that analysis to a future study. Including reionization physics has been shown to significantly boost the number of quiescent galaxies at earlier times (z ∼ 5.5; H. G. Chittenden et al. 2025), which could ultimately become important in resolving tension with future theoretical models.

In the future, the lack of massive quiescent galaxies in state-of-the-art galaxy formation simulations needs to be investigated further to separate two compounding issues. One is the potential overall lack of massive galaxies (either star forming or quenched; e.g., A. Weibel et al. 2024; M. Shuntov et al. 2025), which would point to star formation not being efficient enough or, conversely, feedback being too strong in regulating star formation in the early Universe. The second one is AGN feedback itself, and whether the processes it encompasses (e.g., mechanical, radiative, or energetic feedback) are sufficient to quench massive galaxies in the early Universe. It is clear that this field is nascent, and further observations of massive quiescent galaxies and their stellar mass distribution over larger samples would provide invaluable constraints for galaxy formation models.

A final resolution is empirical: if most apparently quiescent galaxies host heavily dust-obscured star formation, then the apparent tension could disappear. Testing this would require additional observations of apparently quiescent galaxies at 3 < z < 5 in the mid-IR (MIR) or far-IR (FIR). All of the spectroscopic identifications of quiescent galaxies thus far rely on interpreting their rest-frame optical–NIR emission. These inferences cannot yet rule out extreme birth-cloud dust attenuation (A_V ∼ 5) that could hide instantaneous star formation. This scenario has already been discovered in some z < 1 optically selected poststarburst galaxies, which have FIR SFRs ∼ 1–2 dex higher than those inferred from their optical information (D. Baron et al. 2023). As revealed by our PCA (also in O. R. Cooper et al. 2025), the prevalence of galaxies that simultaneously host evolved stellar populations while being dust attenuated suggests this is plausibly a more common scenario at z > 3. If ∼90% of the rest-frame optically selected quiescent galaxies at z > 4 turned out to host star formation embedded in optically thick dust, the 1 dex discrepancy would disappear. Stellar population synthesis modeling of a truly panchromatic sample of quiescent galaxies at z > 3 could lay this uncertainty to rest; the attenuated radiation from instantaneous star formation would inevitably reradiate at MIR and FIR, testable by deeper-than-existing observations with facilities such as the Atacama Large Millimeter/submillimeter Array.