The Formation Histories of Massive and Quiescent Galaxies in the 3 < z < 4.5 Universe

FAST++ is a spectral fitting code written in C++, with similar functionality to IDL FAST code M. Kriek et al. (2009). The SED fitting procedure utilized by us here largely mirrors that of C. Schreiber et al. (2018a) with a few exceptions and upgrades. FAST++ v1.5.0 implements an LSF functionality, where the model spectra can be convolved to a user-defined wavelength dependent σ of the Gaussian LSF. This can be provided using the SPEC_LSF_FILE option. In terms of the stellar population properties, similar to C. Schreiber et al. (2018a), we use G. Bruzual & S. Charlot (2003) high-resolution stellar population models with Padova 1994 stellar tracks (G. Bertelli et al. 1994) with STELIB (J.-F. Le Borgne et al. 2003) and BaSeL v3.1 (P. Westera et al. 2002) spectral libraries. We use a G. Chabrier (2003) IMF and a D. Calzetti et al. (2000) dust law following what was used by C. Schreiber et al. (2018a). In C. Schreiber et al. (2018a), the stellar metallicity of the models was kept at Z_⊙. In this analysis, we are performing spectral fitting utilizing the full rest-frame optical wavelengths. Thus, to take optimal use of available data, we allow the stellar metallicity of the models to vary between the full available grid, which is Z = 0.004, 0.008, 0.02, 0.05 (20%–250% Z_⊙). When transitioning from a fixed-metallicity run to one where metallicity is allowed to vary, the best-fit SFHs can change due to degeneracies between age and metallicity (C. Conroy 2013). Therefore, direct one-to-one comparisons of the SFHs between our results and those of C. Schreiber et al. (2018a) are not possible.

The SFH parameters used in the fitting are also kept the same as C. Schreiber et al. (2018a), which is specifically suited to model quiescent galaxies at z > 3. The model use two epochs to describe the SFHs. The main component includes an exponentially increasing and decreasing SFH, where the two exponents (e-folding times; τ_rise and τ_decl) are kept as free parameters. Additionally, the lookback time that defines the boundary between exponentially increasing and decreasing (t_burst) is also allowed to vary freely. This allows the flexibility in the SFH to capture the slow/fast rising SFHs, constant SFHs, and slow/fast quenching SFHs. Thus, the primary mode of the SFH can be defined as follows:

$$\begin{eqnarray}{\mathrm{SFR}}_{\,\rm{Base}\,}(t)\propto \left\{\begin{array}{cc}{e}^{({t}_{\,\rm{burst}\,}-t)/{\tau }_{\,\rm{rise}\,}} & \,\rm{for}t\gt {t}_{\rm{burst}\,},\\ {e}^{(t-{t}_{\,\rm{burst}\,})/{\tau }_{\,\rm{decl}\,}} & \,\rm{for}t\leqslant {t}_{\rm{burst}\,}\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where SFR_Base(t) is the base SFR at lookback time t. However, this base SFR is unable to capture the most recent star formation episode of a galaxy to where IR and submillimeter observations are sensitive to. In particular, stacked Herschel and ALMA measurements are unable to be recovered by this base SFR parameters. In order to alleviate this, an SFR multiplicative factor, R_SFR, was introduced to the model to be activated close to the time of observation of the galaxy. The time window for which this multiplicative factor was applied (t_free) was allowed to vary freely between the last 10–300 Myr of the galaxies SFH. Thus, the final SFR of the galaxy at lookback time t (SFR(t)) can be mathematically expressed as:

$$\begin{eqnarray}\mathrm{SFR}(t)={\mathrm{SFR}}_{\,\rm{Base}\,}(t)\times \left\{\begin{array}{cc}1 & \,\rm{for}\,\,t\gt {t}_{\,\rm{free}\,},\\ {R}_{\mathrm{SFR}} & \,\rm{for}\,\,t\leqslant {t}_{\,\rm{free}\,}.\end{array}\right.\end{eqnarray} \tag{ 2 }$$

We refer the readers to C. Schreiber et al. (2018a, 2018b) for further detailed justification of this SFR within the context of massive and quiescent galaxies at z > 3.

Table 3 details a summary of the stellar population and SFH parameters used in this analysis. All SFR related priors are varied in logarithmic steps as outlined in the table. Observed photometry is fit simultaneously with the observed spectra. G. Bruzual & S. Charlot (2003) models included with FAST++ do not include emission-line contributions. Thus, for galaxies with strong emission lines (such as Hβ, [O iii]λ5007, Hα, [S ii]λ6716λ6731), we remove regions where the strong emission lines fall as shown by Figures 1 and 2 and broadband photometry contaminated by the strong emissions lines is also removed. For galaxies with evolved stellar populations, the contribution of the nebular continuum to the optical flux is expected to be negligible compared to the stellar continuum (N. Byler et al. 2017). This effect strongly depends on stellar age, with younger stars producing a more prominent nebular continuum. The impact is particularly noticeable in the UV range but diminishes significantly at wavelengths longer than the Balmer break, where the nebular continuum becomes much weaker. The best-fit spectrum is considered as the model with the lowest reduced χ² value, and its values are considered to be the best-fit model parameters. One-thousand Monte Carlo simulations are performed to obtain the 68%, 95%, and 99% confidence intervals of all of the parameters obtained through FAST++ fitting. Both the observed photometry and spectral fluxes are permutated within their respective 1σ errors in each of the Monte Carlo iterations.

Table 3. Stellar Population and SFH Parameters Used in FAST++ Spectral Fitting

Free Parameter	Lower Bound	Upper Bound	Step Size
tburst (Gyr)	0.01	tH(z)	${\mathrm{log}}_{10}(0.05)$
τrise (Gyr)	0.01	3	${\mathrm{log}}_{10}(0.1)$
τdecl (Gyr)	0.01	3	${\mathrm{log}}_{10}(0.1)$
RSFR	10−2	105	${\mathrm{log}}_{10}(0.1)$
tfree (Myr)	10	300	${\mathrm{log}}_{10}(0.5)$
A(V)	0	6	0.1
Free Parameter	Values
Z	0.004, 0.008, 0.02, 0.05
Fixed Parameter	Value
SSP Model	G. Bruzual & S. Charlot (2003)
IMF	G. Chabrier (2003)
Attenuation Curve	D. Calzetti et al. (2000)
z	zspec from JWST/NIRSpec

Download table as:  ASCIITypeset image

The best-fit FAST++ models are shown in Figures 1 and 2. Here, when available, the full range between 1 and 5 μm of the NIRSpec spectrum is used for the fitting. As evident, all galaxies are fit to a very high degree of accuracy with a median reduced χ² of 4.3. 3D-EGS-18996 has the highest reduced χ² of 8.2, which can be partly attributed to the increased flux of the best-fit model ~1 μm. FAST++ is able to model the Balmer/D4000 Å breaks accurately. At these redshifts, these features fall ~2 μm, where the NIRSpec/PRISM resolution is at the minimum and the highest degree of variation. Thus, we found that having an accurate model of the LSF input to FAST++ in the fitting process was crucial to obtain reasonable fits around this wavelength with good overall χ² values.

The SFHs with FAST++ show considerable variety between our galaxies and are shown by Figures 4 and 5. We define the formation time of a galaxy (z_form50) as the redshift at which it formed 50% of its total stellar mass by the time of observation. We define the quenching time (z_quench) as the redshift marking the start of the longest contiguous period before observation during which the SFR fell below 10% of the mean SFR during the galaxy's main formation phase. We consider a galaxy to be a fast quencher if it reaches z_quench within 250 Myr of z_form50; otherwise, we consider it a slow quencher.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3D-UDS-7329 at z = 3.2 is the first galaxy to have reached 50% of its observed stellar mass after the Big Bang, reaching ~6 × 10¹⁰M_⊙ within the first ~600 Myr of the Universe (~1.5 Gyr in lookback time). It is also the first galaxy in our sample to have quenched its star formation by ~1 Gyr of the Universe (~1 Gyr in lookback time). 3D-UDS-35168 also has a similar 50% of observed stellar mass formation time (~700 Myr from the Big Bang); however, it has ~1 dex smaller stellar mass and an extended SFH resulting in it only being quenched ~250 Myr before the time of observation. All of our other sources have only formed ≳1 Gyr after the Big Bang (z < 6).

ZF-COS-10559 is the highest-redshift source in our sample with z = 4.3. It also has the second fastest quenching time at ~1.2 Gyr after the Big Bang. However, given its high redshift, in lookback time it has only quenched within the ~250 Myr before it was observed. In fact, with the exception of 3D-UDS-7329, all of our sources have only quenched within ≲350 Myr of the time of observation. The formation, quenching, and observed times of our sample are further visualized by Figure 6.

Zoom In Zoom Out Reset image size

The best-fit values obtained by FAST++ fitting are provided in Table 4. Since the D. Calzetti et al. (2000) dust law assumes a simplified uniform screen model and lacks the 2175 Å dust bump, we have re-run FAST++ using the M. Kriek & C. Conroy (2013) dust model, finding that the results are statistically consistent with each other.

Table 4. Redshifts and Best-fit FAST++ Parameters of Our Galaxies

Galaxy ID	zspec	M*	Z	Av	SFR10a	zquenchb	zform50c	tSFd	〈SFR〉maine
		${{\rm{log}}}_{10}({M}_{\odot })$		mag	${\mathrm{log}}_{10}({M}_{\odot }\,{\mathrm{yr}}^{-1})$			${\mathrm{log}}_{10}(\mathrm{yr})$	${\mathrm{log}}_{10}({M}_{\odot }\,{\mathrm{yr}}^{-1})$
ZF-COS-10559	$4.30{2}_{4.298}^{4.304}$	$10.41{0}_{10.410}^{10.440}$	$0.05{0}_{0.015}^{0.043}$	$0.30{0}_{0.100}^{0.330}$	$-9.24{0}_{-10.180}^{-5.510}$	$5.01{2}_{4.993}^{5.264}$	$5.23{7}_{5.211}^{6.129}$	$8.16{0}_{7.620}^{8.580}$	$2.28{0}_{1.860}^{2.830}$
ZF-UDS-6496	$3.97{6}_{3.974}^{3.978}$	$10.86{0}_{10.860}^{10.880}$	$0.05{0}_{0.023}^{0.045}$	$0.00{0}_{0.000}^{0.060}$	$-2.51{0}_{-3.910}^{-2.400}$	$4.64{1}_{4.528}^{4.659}$	$5.93{7}_{5.630}^{5.937}$	$7.80{0}_{7.800}^{8.700}$	$3.09{0}_{2.210}^{3.090}$
ZF-UDS-3651	$3.80{7}_{3.806}^{3.809}$	$10.65{0}_{10.640}^{10.660}$	$0.00{8}_{0.004}^{0.007}$	$1.30{0}_{1.200}^{1.330}$	$-1.83{0}_{-2.810}^{-1.170}$	$4.01{2}_{3.984}^{4.012}$	$4.06{9}_{4.033}^{4.118}$	$7.68{0}_{7.200}^{7.690}$	$2.95{0}_{2.950}^{3.440}$
ZF-COS-18842	$3.72{1}_{3.720}^{3.723}$	$10.46{0}_{10.460}^{10.590}$	$0.00{4}_{0.004}^{0.039}$	$0.50{0}_{0.380}^{0.670}$	$-1.90{0}_{-9.880}^{-0.840}$	$4.38{1}_{3.812}^{4.399}$	$4.43{7}_{4.399}^{4.521}$	$7.30{0}_{7.250}^{8.620}$	$3.18{0}_{1.980}^{3.230}$
ZF-COS-20115	$3.71{8}_{3.718}^{3.719}$	$11.08{0}_{11.070}^{11.080}$	$0.05{0}_{0.024}^{0.045}$	$0.50{0}_{0.410}^{0.510}$	$-4.42{0}_{-5.910}^{-2.810}$	$4.09{0}_{4.038}^{4.186}$	$4.75{3}_{4.639}^{4.753}$	$8.45{0}_{8.290}^{8.490}$	$2.65{0}_{2.610}^{2.810}$
ZF-COS-19589	$3.71{7}_{3.714}^{3.721}$	$10.76{0}_{10.740}^{10.760}$	$0.00{4}_{0.004}^{0.016}$	$1.20{0}_{1.060}^{1.190}$	$-9.23{0}_{-15.310}^{-5.680}$	$4.51{5}_{4.235}^{4.561}$	$4.63{7}_{4.341}^{4.720}$	$7.47{0}_{7.400}^{7.680}$	$3.31{0}_{3.090}^{3.370}$
ZF-UDS-4347	$3.70{3}_{3.702}^{3.705}$	$10.45{0}_{10.420}^{10.450}$	$0.02{0}_{0.004}^{0.017}$	$0.40{0}_{0.170}^{0.390}$	$-6.11{0}_{-16.230}^{-4.130}$	$4.04{5}_{4.072}^{4.670}$	$4.45{2}_{4.322}^{4.791}$	$8.36{0}_{7.550}^{8.360}$	$2.11{0}_{2.110}^{2.880}$
3D-EGS-27584	$3.62{4}_{3.622}^{3.626}$	$11.24{0}_{11.230}^{11.250}$	$0.02{0}_{0.010}^{0.019}$	$1.30{0}_{1.210}^{1.350}$	$-10.71{0}_{-10.800}^{-3.790}$	$4.19{9}_{4.115}^{4.215}$	$6.54{6}_{6.204}^{6.675}$	$8.95{0}_{8.920}^{8.960}$	$2.33{0}_{2.310}^{2.360}$
3D-UDS-39102	$3.58{7}_{3.586}^{3.587}$	$10.78{0}_{10.750}^{10.840}$	$0.02{0}_{0.011}^{0.039}$	$1.80{0}_{1.660}^{1.870}$	$0.04{0}_{-1.280}^{0.450}$	$3.64{2}_{3.642}^{3.684}$	$3.74{8}_{3.702}^{3.752}$	$7.74{0}_{7.300}^{7.590}$	$3.01{0}_{3.150}^{3.480}$
ZF-UDS-8197	$3.55{0}_{3.550}^{3.550}$	$10.40{0}_{10.390}^{10.410}$	$0.00{8}_{0.004}^{0.007}$	$1.10{0}_{0.970}^{1.090}$	$-3.92{0}_{-6.900}^{-3.870}$	$3.78{8}_{3.783}^{3.863}$	$3.90{5}_{3.887}^{3.974}$	$7.86{0}_{7.410}^{7.930}$	$2.53{0}_{2.470}^{3.010}$
3D-UDS-35168	$3.52{9}_{3.525}^{3.536}$	$10.18{0}_{10.160}^{10.180}$	$0.00{8}_{0.004}^{0.007}$	$0.00{0}_{0.000}^{0.000}$	$-0.27{0}_{-0.520}^{-0.190}$	$4.00{7}_{3.918}^{4.019}$	$7.24{2}_{7.242}^{7.640}$	$8.15{0}_{8.120}^{8.360}$	$2.07{0}_{1.820}^{2.100}$
ZF-COS-20133	$3.48{0}_{3.479}^{3.481}$	$10.34{0}_{10.340}^{10.350}$	$0.00{4}_{0.004}^{0.004}$	$0.70{0}_{0.610}^{0.690}$	$0.16{0}_{-25.320}^{0.530}$	$4.84{3}_{4.720}^{4.983}$	$5.14{3}_{4.983}^{5.201}$	$7.75{0}_{7.600}^{7.930}$	$2.63{0}_{2.470}^{2.790}$
3D-EGS-31322	$3.43{5}_{3.435}^{3.436}$	$10.74{0}_{10.740}^{10.750}$	$0.05{0}_{0.025}^{0.045}$	$0.30{0}_{0.220}^{0.280}$	$-4.06{0}_{-7.560}^{-4.060}$	$3.86{4}_{3.842}^{3.864}$	$3.93{6}_{3.936}^{3.950}$	$7.46{0}_{7.460}^{7.700}$	$3.29{0}_{3.050}^{3.290}$
3D-EGS-18996	$3.25{9}_{3.258}^{3.261}$	$10.88{0}_{10.870}^{10.880}$	$0.05{0}_{0.023}^{0.045}$	$0.00{0}_{0.000}^{0.000}$	$-4.26{0}_{-4.260}^{-2.800}$	$3.41{0}_{3.403}^{3.410}$	$4.03{3}_{3.970}^{4.055}$	$8.74{0}_{8.710}^{8.760}$	$2.16{0}_{2.140}^{2.180}$
3D-EGS-34322	$3.22{6}_{3.223}^{3.230}$	$10.18{0}_{10.080}^{10.190}$	$0.02{0}_{0.007}^{0.018}$	$1.90{0}_{1.750}^{2.020}$	$1.14{0}_{0.190}^{1.230}$	$3.25{3}_{3.240}^{3.271}$	$3.39{3}_{3.314}^{3.415}$	$7.44{0}_{7.430}^{7.830}$	$2.71{0}_{2.260}^{2.720}$
ZF-UDS-7329	$3.20{7}_{3.206}^{3.209}$	$11.10{0}_{11.060}^{11.120}$	$0.02{0}_{0.010}^{0.018}$	$0.30{0}_{0.140}^{0.270}$	$-2.44{0}_{-20.520}^{-0.980}$	$5.75{8}_{5.219}^{6.390}$	$8.76{0}_{8.086}^{12.318}$	$8.38{0}_{8.340}^{8.710}$	$2.77{0}_{2.440}^{2.780}$
ZF-UDS-7542	$3.20{4}_{3.204}^{3.206}$	$10.66{0}_{10.640}^{10.680}$	$0.05{0}_{0.022}^{0.045}$	$1.40{0}_{0.990}^{1.370}$	$0.69{0}_{-2.370}^{0.850}$	$3.41{4}_{3.314}^{3.424}$	$3.47{2}_{3.472}^{3.699}$	$7.41{0}_{7.410}^{8.780}$	$3.25{0}_{1.920}^{3.250}$
3D-UDS-41232	$3.12{1}_{3.119}^{3.122}$	$11.17{0}_{11.160}^{11.170}$	$0.05{0}_{0.021}^{0.045}$	$0.40{0}_{0.310}^{0.390}$	$-1.18{0}_{-2.150}^{0.100}$	$3.50{2}_{3.502}^{3.629}$	$4.21{4}_{4.119}^{4.214}$	$7.91{0}_{7.850}^{8.160}$	$3.29{0}_{3.080}^{3.350}$
ZF-COS-14907	$2.99{9}_{2.998}^{3.000}$	$10.59{0}_{10.580}^{10.600}$	$0.05{0}_{0.024}^{0.045}$	$0.30{0}_{0.180}^{0.290}$	$-6.57{0}_{-6.490}^{-1.510}$	$3.23{5}_{3.218}^{3.287}$	$3.83{7}_{3.837}^{4.064}$	$8.72{0}_{8.720}^{8.860}$	$1.91{0}_{1.770}^{1.910}$

Notes. 1σ upper and lower bounds derived using Markov Chain Monte Carlo iterations for slinefit spectroscopic redshift measurements and FAST++ inferred values are included with the best-fit parameters.

^aThe SFR of the galaxy computed over a lookback time of 100 Myr. ^bThe redshift at which the galaxy is considered quenched. ^cThe redshift at which the galaxy formed 50% of its total stellar mass at the time of observation. ^dThe length of time where 68% of the total integrated SFR of the galaxy took place. ^eThe average SFR in the time window defined by t_SF.

Download table as:  ASCIITypeset image

Next, we use the Prospector SED fitting code (B. D. Johnson et al. 2021) to infer the SFHs of our galaxies utilizing a flexible approach to parameterize the SFHs. We use a continuity_sfh to parameterize the SFHs of our galaxies. While this is commonly known as a nonparametric SFH, in reality the inference of the SFH is made on fixed time bins in lookback time for each galaxy. The amount of star formation that can happen in each time bin is not constrained; thus, this allows for greater freedom in defining the formation history of galaxies. The SFR in each time bin represents the average star formation happening in that time bin. The SFH for our galaxies is constrained over seven SFR bins, similar to J. Leja et al. (2019a). In lookback time, for each galaxy, the first two bins are fixed between 0–30 Myr and 30–100 Myr, and the final time bin is fixed to be between 85% and 100% of the time of the Universe. The remaining four bins are distributed evenly between the 100 Myr and 85% of the time of the Universe in equally spaced logarithmic time bins. We note that the time bins are slightly different to what was utilized for 3D-UDS-7329 in K. Glazebrook et al. (2024), which was driven by the rapid formation of that source within the first few hundred million years of the Universe.

We allow the stellar metallicity, dust optical depth, and the stellar mass to vary freely. The stellar metallicity is allowed to vary between $-2.5\lt {\mathrm{log}}_{10}(Z/{Z}_{\odot })\lt 0.5$, the dust optical depth following D. Calzetti (2000) dust law is allowed to vary between 0 < τ_dust < 4.0, and the stellar mass is allowed to vary between 7.0 < log₁₀(M_*/M_⊙) < 12.0. The SFR is allowed to vary freely between seven time bins (j) and is parameterized using six log₁₀(SFR) ratios defined as log₁₀(SFR_j/SFR_j+1). We use a G. Chabrier (2003) IMF, MILES stellar library (J. Falcón-Barroso et al. 2011), and MIST stellar isochrones (B. Paxton et al. 2015; J. Choi et al. 2016; A. Dotter 2016). To aid with direct comparisons with FAST++ we turn off nebular line and continuum emission and AGN contribution in Prospector. The effect of nebular contribution is considered in Section 3.6.

Similar to FAST++, we run Prospector using all available photometry (either from JWST/NIRCam or 3DHST/ZFOURGE) and trim the observed spectra at <0.7 μm rest frame. We further mask out emission lines in the observed spectra using the same masks as used with FAST++. The spectral fits agree well with the observations, with all fits showing a mean reduced χ² ~ 3.9 ± 1.7. We find that galaxies with ≳1.5 Z_⊙ tend to have slightly higher χ² values.

Figure 7 shows the stellar mass estimates obtained with FAST++ and Prospector. All masses agree within a ~0.4 dex scatter, which is consistent with the expected scatter between FAST++ and Prospector stellar mass estimates (J. Leja et al. 2019b).

Zoom In Zoom Out Reset image size

The reconstructed SFHs of the galaxies are shown in Figure 4. Driven by the parametric nature of the FAST++ SFHs, the buildup of stellar mass is gradual over time. In comparison, Prospector has more variability in the SFH driven by the increased freedom in its SFH parameterization. While most galaxies exhibit good overall agreement in their SFH shapes between Prospector and FAST++, the parametric nature of FAST++'s SFH modeling results in more narrowly defined star formation burst episodes. Visually, ZF-COS-10559 and ZF-UDS-3651 show the largest deviations in their SFHs between FAST++ and Prospector. ZF-COS-10559, according to Prospector, has an early intense starburst phase that makes ~99.9% of the total formed mass. As we discuss in Section 3.7, ZF-UDS-3651 is a potential AGN with strong emission lines, and ZF-COS-10559 hits the highest possible metallicities in both FAST++ and Prospector fitting. Both of these could add extra complexity to the spectral analysis.

To determine formation time and quenching times as parameterized by the two SED fitting codes, we integrate the cumulative mass of the galaxies in the SFH bins and find the time bin where 50% of the total formed mass at the time of the observation is made. We define this as the formation time window of the galaxy (t_form50). We then find the peak SFR of the galaxy and investigate whether the SFR dips below the necessary 10% threshold of the maximum SFR between the time of the peak SFR to the observed time. If such time exists, we explore whether the galaxies meet this conditions continuously between the peak SFR and the time of observation. If so, we consider the galaxy to be quenched.

In Figure 5, we compare the formation and quenching time of our sample for FAST++ and Prospector SED fitting codes. Given the parametric nature of the SFH, FAST++ reports a best-fit formation and quenching times along with their Monte Carlo uncertainties. Prospector utilizes fixed time bins to infer the variation between the SFR ratios in adjoining time bins; thus, the formation and quenching time are defined in terms of the time bins and are highlighted in blue and orange, respectively. While the formation and quenching time windows agree well between the two codes for most of our sources, there is a general tendency for Prospector to have older t_form50 times (50% of the mass formed at a higher-z) compared to FAST++. Nonparametric SFHs from Prospector have been shown to point toward older SFHs compared to parametric SFHs in the z < 3 Universe (e.g., J. Leja et al. 2020). Our assumption of a flat prior for the continuity_sfh in Prospector differs from the exponentially increasing SFH parameterization assumed by FAST++ at earlier times. Due to the flat prior, Prospector can place more mass formation at earlier times. We further investigate the impact of the assumed SFH prior shape in Section 3.5.

For galaxies that satisfy the quenching criteria in Prospector, the agreement between the FAST++ and Prospector quenching time is better than the formation time. One source, 3D-UDS-39102, is not classified as quiescent based on Prospector results. 3D-UDS-39102 has the lowest S/N in our sample and shows only marginal evidence for a Balmer break based on NIRCam photometry. FAST++ results also suggest very recent quenching for this source (<100 Myr from the time of observation). C. Schreiber et al. (2018a) also found 3D-UDS-39102 to have a lower bound for quenching time to be consistent with 0 Myr with FAST++; however, our joint spectra + photometry fit from FAST++ suggest a 1σ lower bound of ~​​​​​​30 Myr for quenching time. 3D-EGS-34322 satisfies the quenching criteria only in the most recent SFH bin in Prospector, which is defined between 0 and 30 Myr. This aligns well with the FAST++ quenching time of approximately 20 Myr.

Our results demonstrate that the parametric form of the SFH optimized for z > 3 massive galaxies presented by C. Schreiber et al. (2018a) largely provides consistent results to the SFHs parameterized using Prospector with a more flexible approach. However, there are other notable assumptions routinely made in SED fitting such as the choice of input stellar population models, how emission lines are considered, and how effects of an AGN are considered in SED fitting. For example, the MIST stellar library used with Prospector includes effects of stellar rotation while the BC03 stellar population models we used for FAST++ do not consider these effects. Thus, it is imperative to understand the role different SED fitting assumptions such as the input stellar libraries, emission line contributions, and AGN contributions play in determining the formation histories of our massive quiescent galaxies. Given Prospector by design affords flexibility to vary these parameters, we opted to use Prospector for the next steps of our analysis and only included galaxies that have full rest-frame optical coverage in our observed spectra. A brief summary of the Prospector assumptions that we utilized is presented in Table 5.

Table 5. Summary of SED Fitting Assumptions Explored Using Prospector

Free Parameter	Range
SFH bins (lookback time in gigayears)	0–30, 30–100, bin3–6, 0.85tuniv- tuniv
Dust optical depth (D. Calzetti 2000)	0 < τdust < 4.0
Stellar mass	$7.0\lt {\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })\lt 12.0$
Ionization parametera	−4.0 < U < −1.0
AGN bolometric luminosity fractiona	10−5 < fagn < 3.0
AGN optical depth for individual cloudsa (M. Nenkova et al. 2008)	5 < τagn < 150
Marginalize emission linesa, b (B. D. Johnson et al. 2021)	100 < σ < 20,000 km s−1
MILES stellar library (J. Falcón-Barroso et al. 2011) & MIST stellar isochrones (B. Paxton et al. 2015; J. Choi et al. 2016; A. Dotter 2016)
Stellar Metallicity	−2.5 < log10(Z/Z⊙) < 0.5, Z⊙ = 0.0142
C3K stellar library (C. Conroy et al. 2009) & MIST stellar isochrones
Stellar metallicity	−2.5 < log10(Z/Z⊙) < 0.5, Z⊙ = 0.0142
BPASSv2.0 stellar population models (J. J. Eldridge et al. 2017)
Stellar metallicity	−2.3 < log10(Z/Z⊙) < 0.3, Z⊙ = 0.0200
Fixed parameter	Value
IMF	G. Chabrier (2003)
SFH prior	continuity_sfh with a flat priorc
Spectral smoothing (smooth_type)	LSF
emcee sampling (D. Foreman-Mackey et al. 2013)	nwalkers=1024, nburn=[16, 32, 64], niter=2056

Notes.

^aRuns are also computed completely turning off this parameter. ^bOnly the following lines are marginalized over: Hβ, [O iii]λ4959λ5007, [N ii]λ6548λ6584, Hα, [S ii]λ6716λ6731. ^cResults are also compared with SFR(t) ∝ e^t/τ prior with τ = t_univ/4, where t is time in the age of the Universe, and t_univ is the age of the Universe corresponding to the spectroscopic redshift of the galaxy (J. Leja et al. 2019a).

Download table as:  ASCIITypeset image

We investigate three different stellar population models that are inbuilt with FSPS (C. Conroy et al. 2009). First, we use the default MILES stellar library (J. Falcón-Barroso et al. 2011) with MIST stellar isochrones (B. Paxton et al. 2015; J. Choi et al. 2016; A. Dotter 2016) to generate the model spectra used by Prospector to make inferences on the observed spectra and photometry of our sources. MILES library has a higher resolution of R ~ 3000 in the optical and a lower resolution at λ > 7500 Å. Second, we use the C3K stellar library used in FSPS (C. Conroy et al. 2009). Here the spectral templates are sampled at a R ~ 3000 between 2750 < λ < 9100. At z ~ 3.0, this cutoff translates to 3.64 μm in the observed frame. Finally, we also consider BPASSv2 stellar population models (J. J. Eldridge et al. 2017). BPASS spectral resolution is comparable to C3K, but has effects of binary stellar evolution incorporated in stellar evolution. BPASS spectra are sampled at 1 Å uniform sampling, but underlying stellar atmospheric models will have a lower resolution, possibly limited by the resolution of the BaSeLv.31 stellar library.

Fitting is performed with similar parameters as explained in Section 3.2 using the three spectral libraries/stellar population models. SFHs reconstructed using the MILES stellar library (J. Falcón-Barroso et al. 2011) are shown in blue, the C3K stellar library (C. Conroy et al. 2009) in green, and the BPASS stellar population models (J. J. Eldridge et al. 2017) in purple. The residuals of the best fits are shown in Appendix Figure A1. Observed spectra are trimmed at rest frame <0.7 μm. Prospector recovered maximum a posteriori SFHs are shown in Figure 8. The reduced χ² of spectral fits between the three models agree well within 1σ of each other. All of our sources show largely consistent formation histories between the three models. ZF-UDS-4347 is a notable exception, where the C3K models prefer a very early formation while the other models prefer a later growth.

Zoom In Zoom Out Reset image size

We further explore the variation in results between these three models using the definitions for galaxy formation and quenching timescales as discussed in Section 3.2. We show the reconstructed formation and quenching timescales for our galaxies in Figure 9. The formation and quenching times between the three models largely agree with each other. The notable exception is ZF-UDS-4347, where the C3K models prefer an earlier formation time compared to the other models. 3D-UDS-39102 is considered not quenched by all three models. Similarly, the quenching time bin for 3D-EGS-34322 is also in the latest time bin for all three models. There are no instances where the choice of model determines whether a galaxy is classified as quenched or not.

Zoom In Zoom Out Reset image size

Galaxy evolution is a statistical study. Galaxy evolution models by design are tuned to reproduce a diverse range of galaxy scaling relations that describe the average properties of the observable Universe (R. A. Crain & F. van de Voort 2023). Therefore, it is plausible that the models find it challenging to reproduce the properties of individual sources that are outliers from the general galaxy population. In Figure 8, we have shown that the input stellar population models used in the fitting can have some effect on SFH recovery in an individual galaxy basis. In the context of our analysis, we argue that the average formation histories of massive galaxies are more important than individual reconstructions of their formation histories.

To investigate this further, we explore if the three different stellar libraries/stellar population models used in our analysis provide consistent results to the average formation history of our sources. We resample the SFR of each galaxy in 1 Myr increments up to 2200 Myr and compute the average SFR at each time step for our sample, weighting by the number of galaxies in each time step. The 1σ scatter at each time step, calculated using the normalized absolute deviation, is considered as the uncertainty in the SFH. In Figure 10 we show the evolution of the average SFH for our galaxies as a function of cosmic time. The average SFR evolution across the three models is in good agreement. All models show a consistent increase in SFR for ages greater than 1 Gyr, peaking around 1.5 Gyr. This peak is followed by a decline in the average SFHs due to a reduced number of galaxies beyond 1.5 Gyr and the fact that our sample consists of massive quiescent galaxies with limited SFRs in the last few hundred megayears. All models show an enhancement of SFR in the first ~250 Myr, with the C3K model exhibiting the largest increase.

Zoom In Zoom Out Reset image size

From Figure 10, it is clear that regardless of the spectral library/stellar population model used with Prospector, the average SFHs are statistically consistent with each other. The major deviation is only observed in the first ~250 Myr of the Universe. The NIRSpec prism largely covers the rest-frame optical regime. With continued star formation, features from younger stellar populations are likely to dominate the observed spectrum. This, combined with the lower spectral resolution, means that spectral fitting will have less sensitivity to F and G star features that would indicate the nature of the oldest stellar populations in a galaxy. Therefore, a larger uncertainty is expected in this first ~250 Myr window, which may translate to deviations between the different models.

To investigate the effect of utilizing the rest-frame NIR in spectral fitting, we use the Prospector C3K stellar library to fit the spectra of the galaxies, covering the full 1–5 μm spectral range. We utilize C3K for the comparison here because its spectral resolution is higher at longer wavelengths compared to the other models. Specifically, at wavelengths >9100 Å, the critically sampled spectral resolution of the C3K models remains at R ~ 500,¹⁵ which is comparable to the resolution of the observed data (T. Nanayakkara et al. 2024).

The best-fit Prospector models are shown in Figures 1 and 2. The reconstructed average SFHs of the galaxies are shown in Figure 11. The overall shape of the reconstructed SFHs shows good agreement between the λ_rest < 0.7 μm fits and the full λ_observed = 1–5 μm fits. The largest deviation is observed in the earliest time bin, which parameterizes the first ~250 Myr of the Universe. When the full spectral range is utilized in the fitting, the enhancement of the SFR at the earliest times decreases, bringing the average SFR of C3K into close agreement with the SFRs observed for the MILES and BPASS models, as shown by Figure 10.

Zoom In Zoom Out Reset image size

The nonparametric SFH spectral fitting of our analysis performed with Prospector, uses a continuity_sfh model with a flat prior. We utilize seven fixed time bins, as detailed in Table 5. We use smaller time steps for the most recent bins to achieve tighter constraints on recent SFH episodes, while employing logarithmically spaced bins for intermediate epochs to efficiently cover the Universe's age with sufficient resolution. The variability of the SFRs between the bins is parameterized using the logarithmic ratio of the SFR in adjacent bins, defined as:

$$\begin{eqnarray}&&\mathrm{log}\,\rm{SFR\_Ratio}\,={\mathrm{log}}_{10}(\mathrm{SFR}(i)/\mathrm{SFR}(i+1))\end{eqnarray} \tag{ 3 }$$

where i represents the closest time bin in lookback time. A constant prior with a mean of 0, a standard deviation of σ = 0.3, and 2 degrees of freedom, following a Student's T-distribution, is used to sample the posteriors in the Markov Chain Monte Carlo. This resembles a constant SFH prior.

In this Section, we test the role of the assumed SFH prior used in SED fitting. First, motivated by the average accretion rate of halo mass over cosmic time (A. Dekel et al. 2013; C. Turner et al. 2025), we apply an exponentially increasing SFH prior for the $\mathrm{log}\,\rm{SFR\_Ratio}\,$ to explore how it affects the average SFHs of our observed sample. Thus, we modify the constant mean of 0 prior in each of the six $\mathrm{log}\,\rm{SFR\_Ratio}\,$ bins following $\mathrm{SFR}({t}_{i})={e}^{{t}_{i}/\tau }$, with τ = t_univ/4, where t_i is the mean time in the age of the Universe, and t_univ is the age of the Universe corresponding to the spectroscopic redshift of the galaxy (J. Leja et al. 2019a). All other parameters are unaltered. We use the Prospector C3K stellar library to fit the spectra of the galaxies utilizing the updated prior covering the full 1–5 μm spectral range.

The reconstructed average and individual SFHs of the galaxies with flat and exponentially increasing priors are shown by Figures 11 and 12. The overall shapes of the reconstructed SFHs show good agreement between the two priors. Both SFHs exhibit a slight dip in the mean SFR around ~1.25–1.50 Gyr, with the dip occurring at slightly earlier times in the fits using the exponentially increasing prior. The most significant shift between the two SFHs is in the first 250 Myr. The constant SFH prior shows a slight increase in the mean SFR at <250 Myr, while the exponentially increasing SFR shows a slight decrease in the SFR at a similar amplitude. This shift can be attributed to galaxies such as 3D-EGS-18996 and ZF-UDS-7329, where the mass formed in the oldest time bin (the highest-z bin) is reduced by the introduction of the exponentially increasing prior.

Zoom In Zoom Out Reset image size

Similarly, we also investigate the effect an exponentially decreasing prior, $\mathrm{SFR}({t}_{i})={e}^{-{t}_{i}/\tau }$, has on the average SFHs of the observed galaxies. As shown by Figure 11, the average evolution of the SFR > 250 Myr is largely consistent with the flat and exponentially increasing prior. The peak SFRs in all three models are similar, where a gradual incline of the SFRs can be observed from >1 Gyr onward. The major discrepancy is at the first ~250 Myr of the Universe, where the exponentially decreasing SFR shows the highest increase in SFR at the earliest times. From Figure 12, it is evident that an increasing fraction of galaxies tend to prefer more mass formation in the earliest time bin when the exponentially decreasing prior is introduced in Prospector.

Overall, the three assumed SFH priors in Prospector agree well with the FAST++ parametric form at later times. Visually, the exponentially increasing SFH prior shows the closest results to FAST++, which is expected given that FAST++ also utilizes an exponentially increasing parametric form for the SFH. Based on our tests, it is evident that the prior assumed for the posteriori sampling has an effect at older times, where the diagnostic power of the age sensitive features may be limited. This can have an effect on how the earliest stages of the SFH of galaxies are probed with full spectral fitting techniques. At later ages where the observed spectral features show diagnostic sensitivity, consistent results to the SFHs can be obtained independent of the assumed prior.

In both FAST++ and Prospector spectral fitting so far, we have masked out emission lines in the observed galaxies and used models that do not account for contributions from nebular continuum and nebular emission lines. In this Section, we investigate whether allowing nebular contributions as a free parameter in the spectral fitting would significantly affect the reconstruction of the SFHs.

The Flexible Stellar Population Synthesis code (FSPS; C. Conroy et al. 2009), utilized by Prospector, computes emission lines using the cloudy photoionization code (G. J. Ferland et al. 2017), which are stored in a precomputed grid (N. Byler et al. 2017). The age and gas-phase metallicity of the input spectra are matched to generate nebular continuum and emission lines for the FSPS stellar population model. This combined spectrum is then used for inference with the observed data. The emission lines are powered solely by the ionizing photons from the stellar population models, with no AGN contribution included in the ionizing spectrum.

To account for possible contributions from nebular emission lines, especially given the high-EW limits observed in the NIRSpec prism mode, we marginalize over the cloudy input grid with variable ionization parameters. We remove the emission line masks applied to the observed spectra and allow for emission line contributions in Prospector by letting the gas ionization parameter vary freely from −4.0 < U < −1.0. The gas-phase metallicity is fixed to match the stellar metallicity. Each line is fitted with a single-component Gaussian and convolved with the instrument's LSF.

Figure 13 shows a comparison of the formation and quenching time for the galaxies when Prospector is run with and without nebular contribution as a free parameter. We note that when nebular contributions are not considered, the emission lines in the observed spectra are masked. For the majority of galaxies, both runs produce similar formation and quenching times. The largest deviation in formation time is observed for 3D-EGS-34322, where the formation time window is approximately 500 Myr earlier than when nebular contributions are considered. A similar offset in quenching time is seen for ZF-UDS-3651. In both runs, 3D-UDS-39102 is not considered quenched. When emission line contributions are included, 3D-UDS-8197 is also classified as not quenched by Prospector. The Prospector SFH reconstruction of 3D-UDS-8197 indicates that the galaxy was classified as quiescent during a period before the observation but has since experienced an increase in SFH, classifying it as star-forming at the time of observation. The rest-frame UVJ colors, however, classify this galaxy as quiescent, which is discussed further in Section 4.2. Additionally, this galaxy exhibits strong and broad emission lines likely driven by AGN activity. Given that the cloudy photoionization models used in Prospector are solely driven by star formation, they tend to elevate the SFR to match the observed spectrum, thereby failing to account for possible AGN contributions to the emission lines. The limitations of this approach are further discussed and addressed in Section 3.7.

Zoom In Zoom Out Reset image size

In this Section, we explore the role of adding AGN contribution as a free parameter to SED fitting for SFH recovery. There are two free parameters that are introduced to the Prospector fitting to account for AGN contributions. As summarized in Table 5, the first is the AGN bolometric luminosity fraction, which is defined as a fraction of the stellar bolometric luminosity. This is allowed to vary from 10⁻⁵ < f_agn < 3.0. The second parameter related to AGNs is the optical depth for individual clouds, as parameterized by M. Nenkova et al. (2008). This is computed in the V band and is allowed to vary from 5 < τ_agn < 150. The contribution of AGNs to the emission lines is not accounted for in the Prospector fitting. Consequently, the emission lines are not physically motivated by an underlying stellar or AGN ionizing spectrum, and thus the SFH is associated solely with the stellar continuum. Effects of AGNs will be addressed further in a forthcoming paper by M. Martínez-Marín et al. (2025, in preparation).

Based on Figure 14, for most galaxies, we find no statistically significant deviations in the reconstructed SFHs when AGN contributions are included. However, considering the overall shape of the SFH, both ZF-COS-10559 and ZF-COS-19589 show notable deviations between the two runs. In ZF-COS-10559, once AGN contributions are considered, the SFH shows a gradual incline, reducing the contribution of an early burst to the final stellar mass. Conversely, in ZF-COS-19589, the AGN contributions suggest that more mass was formed at earlier times.

Zoom In Zoom Out Reset image size

There are six galaxies in our sample with either significant detections of Hα or [O iii]λ5007 emission lines. We define significant detection as a line with S/N > 5 and an observed line EW of  > 10 Å, and we show these galaxies in Figure 15. Visual selection of the spectra also confirms the lines of these sources to be clearly visible over the local continuum level. In ZF-UDS-3651 and ZF-UDS-8197, both Hα and [O iii]λ5007 satisfy the significant detection criteria; in others, only one of the two lines satisfies this selection. Visual inspection of ZF-UDS-7542 might suggest Hα to also be significantly detected; however, it is observed at an EW ~ 6 Å, and thus is lower than our EW cut. With a nominal R ~ 200 in PRISM mode observations, our velocity resolution is ≳1500 km s⁻¹ in the ~K band. Therefore, we are unable to provide tight constraints to the observed line widths of these sources, except for 3D-EGS-18996. Both [O iii]λ5007 and Hα are broad for this galaxy with slinefit providing a best-fit line width of ~6000 km s⁻¹.

Zoom In Zoom Out Reset image size

The ionizing spectrum of the Prospector cloudy grids are limited to stars; hence, emission lines lack contribution from AGN sources. While there are only minimal changes to the recovered SFHs when AGN effects are added, it is possible that this limitation of the input cloudy photoionization grids has an effect on the Prospector results. Given the prominence of the emission lines of the six galaxies highlighted in Figure 15, it is possible for them to have a contribution from an AGN. While we cannot completely rule in favor of strong AGN contribution for these sources, it is necessary to determine the limitations of using cloudy grids that only have contributions from stars to our recovered SFHs.

We explore this by removing the emission line masks in the observed spectra and allowing freedom in Prospector to fit emission lines independently of the input ionization conditions. In this setting, emission lines are considered purely Gaussian in nature and are marginalized over line width prior to obtain the maximum a posteriori solution. We utilize the nebular_marginalization template in Prospector and select the following strong optical emissions to marginalize over: Hβ, [O iii]λ4959λ5007, [N ii]λ6548λ6584, and Hα, [S ii]λ6716λ6731. Given the prism resolution, Hα and [N ii]λ6548λ6584 are not resolved.

As shown in Figure 14, the recovered SFH reconstructions are largely statistically consistent between the two Prospector runs that consider AGN effects. We also find the f_agn and τ_agn parameters to be statistically consistent between the two runs. When lines are marginalized, f_agn varies between 0.25 and 0.93, with ZF-UDS-8197 showing highest f_agn value. Similarly, τ_agn varies between 5.9 and 7.9, with 3D-EGS-18996 showing the largest value. However, we note that constraining both the ionizing spectrum and the overall SED shape together with star-forming and AGN contributions is necessary and would provide tighter constraints to the nature of the AGN in galaxies such as ours. We revisit this in our forthcoming paper M. Martínez-Marín et al. (2025, in preparation) and envision that novel advancements in machine-learning-assisted spectral fitting techniques would be able to address this in future (e.g., Y. Li et al. 2024).

So far we have investigated the role of the form and prior of the assumed SFH parameterization, the input stellar library/stellar population model, the wavelength window of the observed spectrum used for spectral fitting, nebular emission, and AGNs in reconstructing the SFHs of our massive quiescent galaxy sample. The results indicate that the final best-fit or maximum a posteriori SFH solution in general has good agreement between these variations. The average SFH of the sample is also shown to be statistically consistent. However, the main limitation here is that we have no ground truth to establish the accuracy of the recovered SFHs.

To assess this limitation, we utilize mock observations to investigate whether the input SFHs can be recovered accurately by Prospector. If Prospector works as intended, the input SFHs of our models should be recovered through the fitting. This is because Prospector utilizes a comprehensive Bayesian framework that effectively leverages the available wavelength coverage and S/Ns to constrain the SFH parameters. By accurately modeling the observed spectral energy distributions and accounting for uncertainties, Prospector is expected to reproduce the input SFHs when the underlying assumptions and observational conditions are met. In order to use realistic SFHs to perform mock observations, we utilize the IllustrisTNG simulations (A. Pillepich et al. 2018a; V. Springel et al. 2018) to select 283 massive galaxies (M > 10¹¹M_⊙) at z = 3 from the TNG300 suite that are expected to be quiescent based on the definition of specific star formation rate (sSFR) < 1.5 × 10⁻¹⁰ yr⁻¹ from C. Schreiber et al. (2018a). We note that TNG300 simulations show good agreement for the number density of massive quiescent galaxies at z ~ 3–4 (F. Valentino et al. 2023).

Once the galaxies are selected from TNG300, the SFHs are computed from the mass-weighted distribution of stellar formation times of all star particles bound to the subhalo, as was done in H. G. Chittenden & R. Tojeiro (2023). This gives the SFHs a greater time resolution than the snapshot data: 1.4 Myr per age bin, compared with a median snapshot time difference of 74 Myr, up to z = 3. Furthermore, unlike the stellar mass obtained from the merger tree, this field is derived from the stellar ages of all particles of the subhalo. Therefore, we can compute an intrinsic SED from a fine-grained composite stellar population. Additionally, it effectively captures the star formation from all progenitors, not solely the main progenitor branch, therefore tracing the stellar mass contribution of all merger events throughout the galaxy's history. Finally, this results in realistic SFHs for 283 sources.

Mock spectra for the 283 TNG300 selected sources are computed using the FSPS code. The MILES stellar library with MIST stellar isochrones and P. Kroupa (2001) IMF is used for this purpose. Metallicity is fixed at 50% Z_⊙, and a τ_dust = 0.5 is applied to the spectra following a D. Calzetti (2000) dust law. The evolution of the SFR with cosmic time for each galaxy is input into the stellar population code using the set_tabular_sfh option. Spectra are generated at the rest frame and are redshifted to their observed redshift. A physical velocity of 300 km s⁻¹ is applied to the spectra following results from J. Esdaile et al. (2021a). We convert the NIRSpec PRISM instrument LSF to velocity and assert that the intrinsic resolution of the SSPs is lower than that. We then apply the LSF velocity smoothing to the spectra after subtracting the SSP spectral library resolution. Noise is added to the spectra following a random normal distribution such that there is a median S/N of 50, which is similar to our observations. The photometry for the mock spectrum is computed in the following JWST NIRCam bands using the sedpy¹⁶ package: F115W, F150W, F200W, F277W, F356W, F410M, and F444W. For simplicity, we don't add emission lines to the mock spectra nor use it as a free parameter in Prospector fitting.

The mock spectrophotometry is fit using Prospector with similar parameters as outlined in Section 3.2. The recovered average SFHs are shown in Figure 16. The Prospector recoveries show good agreement with the average SFH of the input galaxies, with the growth of the peak of the SFRs t = ~1.5 Gyr and the subsequent quenching phase matched well. The mock tests demonstrate that the earliest time bin defined around the first 300 Myr in the Universe is not well constrained by the SFH recovery. As discussed in Section 3.5, we attribute this offset at the earliest bin to the use of a flat constant SFH prior with Prospector. However, our tests with realistic SFHs from TNG300 massive quiescent sources demonstrates that on average, current SED fitting models are able to recover the SFHs in the z > 3 Universe with good sensitivity to stellar ages, even with low-resolution NIRSpec PRISM data.

Zoom In Zoom Out Reset image size

In order to compare the average formation histories of massive quiescent galaxies in our observed sample with the predictions from cosmological simulations, we utilize four cosmological suites as outlined below. To assist in comparison between the simulations, we use a constant definition of M_* > 3 × 10¹⁰ M_⊙ (mean mass of the observed sample) and sSFR < 1.0 × 10⁻¹⁰ yr⁻¹ as the selection for massive quiescent galaxies.

First, we select galaxies from the TNG300 cosmological suite (A. Pillepich et al. 2018b). TNG300 has a size of 302.6³ Mpc³, and we find 556 sources that are classified as massive quiescent based on our definition at z ~ 3. The SFHs of these galaxies are computed as detailed in Section 3.8. In Figure 17, we show the comparison between the average SFHs of our massive quiescent sample (see Section 4.2) with the average SFHs of the TNG300 massive quiescent sample. Regardless of the stellar population models or the utilized SFH prior, the buildup of the average SFH of the observed sample largely follows what is predicted by TNG300. The significant deviation is during the first ~250 Myr of the Universe. As we noted in Section 3.3, this initial time window also shows the largest dependency for the input stellar library/stellar population models. As discussed in Section 3.5, utilizing an exponentially increasing SFH prior over a constant SFH prior, the average SFHs at the earliest time bins can be reduced. However, even with the exponentially increasing SFH prior, the offsets between the SFHs of the models and the data are not statistically consistent with each other.

Zoom In Zoom Out Reset image size

Utilizing the same IllustrisTNG cosmological suite, we use the SFHs of massive quiescent galaxies from the TNG100 simulations (A. Pillepich et al. 2018b) to compare with the observations. Compared to TNG300, TNG100 has a smaller size of 110.7³ Mpc³ but has ~1 dex higher particle resolution size. Given the limited volume of the TNG 100 simulations, we only find 43 galaxies that satisfy our massive quiescent definition at z ~ 3. We show the average SFHs of these galaxies in Figure 17. While both TNG100 and TNG300 exhibit largely similar mass growth, the SFRs in TNG100 show a slightly earlier rise and fall. Additionally, TNG100 shows a rapid increase in SFR at ~2 Gyr before its rapid decline to z = 3.0.

The lower particle size of TNG100 compared to TNG300 allows for finer resolution in feedback processes and halo mass growth mechanisms (D. Nelson et al. 2018; V. Springel et al. 2018; H. G. Chittenden & R. Tojeiro 2023). The higher particle resolution in TNG100 enables more detailed mapping of the gravitational potential of individual stellar and dark matter particles. As a result, TNG100 produces deeper (higher amplitude) potential wells compared to TNG300. Since star formation is coupled with the simultaneous cooling of gas and an increase in pressure, TNG100 achieves more efficient star formation, leading to the slightly earlier bias in mass growth relative to TNG300. Additionally, the finer resolution of feedback mechanisms in TNG100 allows for more effective regulation of star formation.

Next, we use SHARKv2.0 (C. D. P. Lagos et al. 2024) semianalytical models to compare the formation histories of massive quiescent galaxies at z ~ 3. SHARK has a volume of 464³ cMpc³, and we find 502 galaxies that satisfy our massive quiescent galaxy criteria at z ~ 3. For each galaxy, we investigate the SFR as a function of cosmic time and compute the mean SFH of the sample. We show this in Figure 17. While SHARK galaxies show a marginal increase in the average SFHs at earlier times, similar to the TNG suite, at >250 Myr the shape of mass buildup of SHARK galaxies largely resembles our observed population.

In Figure 17 we further show the average formation histories of z ~ 3 massive quiescent galaxies from the Magneticum Pathfinder hydrodynamical cosmological simulations (R.-S. Remus & L. C. Kimmig 2023). These sources are selected from the snapshot 6 (corresponding to z = 3.4) of Magneticum Box 3 with a boxsize of 128 cMpc h⁻¹. Note that 184 central galaxies satisfy our massive quiescent galaxy criteria. The average mass buildup of the Magneticum sample is similar to the other simulations. However, the Magneticum sample reaches the peak SFRs at earlier times compared to the other simulations. Thus, it is evident that Magneticum has quenching mechanisms that are, in general, switched on at earlier times compared to those of the other simulations. This results in providing a higher number density of massive quiescent sources at z ~ 4 (L. C. Kimmig et al. 2025). We attribute the early decline in the SFR in the Magneticum simulations to the more flexible implementation of feedback mechanisms.

In Magneticum, quenching is attributed to both AGN and star formation feedback, allowing galaxies in underdense regions to expel gas more efficiently and quench star formation (L. C. Kimmig et al. 2025; R.-S. Remus & L. C. Kimmig 2023). Since quenching in Magneticum only requires a rapid isotropic gas inflow to trigger a starburst, a rapid mass buildup may immediately precede quenching. As Magneticum differentiates cold and hot accretion onto the black hole (L. K. Steinborn et al. 2015), this allows for continuous accretion as long as cold dense gas falls in isotropically, even while a hot outflow is starting to be driven by the AGN. This leads to earlier black hole growth and stronger feedback (see L. K. Steinborn et al. 2015). Thus, while AGN feedback plays a significant role in quenching, there is less of a requirement for the black hole to reach a characteristic stellar mass to drive kinetic mode quenching.

The quenching of star formation in TNG300 galaxies at z ~ 4 is attributed to the kinetic mode of the AGN being triggered when the central black hole reaches a critical stellar mass (A. I. Hartley et al. 2023). Once the kinetic mode energy exceeds 1% of the thermal mode energy in the AGN, the kinetic mode AGN expels the remaining star-forming gas in the galaxies, leading to quenching (R. Weinberger et al. 2018). Therefore, in TNG300, black holes must exceed a certain mass threshold before triggering quenching mechanisms, which results in a delayed onset of quenching compared to Magneticum.

Most of the quenching in SHARK happens through the AGN jet-mode feedback, which can start acting as long as jets are produced and there is a hot halo against which to produce work. The latter condition makes it hard for SHARK to start quenching a significant amount of massive galaxies at earlier times, producing the later quenching of galaxies compared to Magneticum.

The normalization difference between the peak SFRs of different simulations is due to the stellar mass distribution of the different galaxies. The mean mass of the TNG suite is ~3×–4× higher compared to SHARK and Magneticum. The observed mean mass of our quiescent galaxies agrees well with SHARK and Magneticum. Thus, the absolute normalization in SFR for our observations is similar to that of SHARK and Magneticum, while TNG100 and TNG300 have a marginally higher normalization, as is evident from Figure 17. This is further shown by Appendix Figure A2.

Further analysis of the formation and quenching mechanisms across a range of hydrodynamical and semianalytical models is presented in C. d. P. Lagos et al. (2025).

Our analysis utilized novel JWST NIRSpec 1–5 μm spectroscopy of ground-based K-selected z > 3 massive quiescent galaxy candidates to explore their formation mechanisms. We showed consistent SFHs between different SED fitting codes, SFH parameterizations, input stellar population models, and other ISM and AGN related properties. Our sample showed a variety of formation and quenching timescales, ranging from the first billion years of the Universe to few megayears before the time of observation. In this Section, we first briefly discuss the advancements made in the z > 3 Universe from previous NIR facilities and then explore the advancements we are now able to make with JWST.

In the pre-JWST era, deep observations from ground NIR facilities provided a first look into the properties of z > 3 massive quiescent galaxies. In K. Glazebrook et al. (2017), we presented the first deep Keck/MOSFIRE K band spectroscopy of a z = 3.7 massive quiescent galaxy, which was photometrically identified by NIR medium band imaging data from the ZFOURGE survey (C. M. S. Straatman et al. 2016). Spectral fitting inferred the galaxy to be ~500–1000 Myr old and to have have assembled ​​​​​​~10¹¹M_⊙within a ≲250 Myr window. This implied a very high star formation efficiency for galaxies at z > 5, which was not observed in direct UV observations of sources at these early times (e.g., R. Smit et al. 2012; B. Salmon et al. 2016). A similar view of the early formation timescales was also established by the study of massive strong emission-line-dominated galaxies at z > 3 (e.g., Z. C. Marsan et al. 2017).

In C. Schreiber et al. (2018a), we preformed a mass complete analysis of the z > 3 massive quiescent candidates utilizing photometric selections from ZFOURGE (C. M. S. Straatman et al. 2016) and 3DHST (R. E. Skelton et al. 2014) surveys. Keck/MOSFIRE K- and H-band spectroscopy of 24 candidates was presented by C. Schreiber et al. (2018a). Twelve sources were spectroscopically confirmed, out of which two were found to be z ~ 2 interlopers. The remaining sources, which primarily constituted the fainter end of the K-band selection eluded confirmation (also see B. Forrest et al. 2020b), even with deep ~10 hr of Keck/MOSFIRE observations. A reconstruction of the SFHs of the 22 candidates showed a variety of formation and quenching timescales, but stringent limits were not possible due to the limited S/N and the limited spectral coverage of the ground-based data.

Utilizing deep KECK/MOSFIRE and VLT/XSHOOTER observations, many surveys of z > 3 massive quiescent galaxies were conducted by several teams. This result transformed our view of rapid formation of massive sources in the early Universe. The MAGAZ3NE survey presented spectroscopic confirmation for 16 photometrically selected >10¹¹M_⊙ galaxies at z > 3 selected from the UltraVISTA survey (B. Forrest et al. 2020b), out of which ~5 were expected to be quiescent. The ages of the massive quiescent candidates were constrained using the D4000 Å and Hδ spectral features, which suggested that most of the stellar mass of these galaxies formed <1 Gyr from the time of observations. Tighter constraints were obtained with multiband SED fitting, which constrained the formation window for MAGAZ3NE quiescent sources to be between 4 < z < 5 with intense star formation up to ~3000M_⊙ yr⁻¹ followed by rapid quenching. F. Valentino et al. (2020) presented spectroscopy of three ~10¹¹M_⊙ quiescent galaxies at z ~ 3.7–4.0. They found the galaxies to have likely formed most of their stellar masses by z ~ 5 and subsequently quenched at z ≲ 5 after experiencing peak SFRs of ~3500M_⊙ yr⁻¹.

A. C. Carnall et al. (2020) presented a sample of 10 robust z > 3 photometrically selected massive quiescent candidates selected from the CANDELS UDS and GOODS-S fields (N. A. Grogin et al. 2011; A. M. Koekemoer et al. 2011). Spectroscopy for the sources was limited to the rest-UV, where Lyα emission was confirmed for two candidates, while another showed evidence for a clear Lyman break (the latter source was spectroscopically confirmed to be a z = 4.7 massive quiescent galaxy; A. C. Carnall et al. 2023b). The reconstructed SFHs for the robust sources suggested that the galaxies assembled most of their stellar masses within the first billion years of the Universe.

J. Antwi-Danso et al. (2025) presented spectroscopy of three >10¹⁰M_⊙ quiescent galaxy candidates at z = 3.3–4.7, selected using split-K band imaging from the FENIKS survey (J. Esdaile et al. 2021b). Two of the sources are expected to have formed at z ~ 4 and undergone subsequent rapid quenching. One of the other galaxies (albeit with very poor spectral quality) suggests a formation time window at z ~ 9. This is similar to what was observed for ZF-UDS-7329 (also see A. C. Carnall et al. 2024; K. Glazebrook et al. 2024; T. Nanayakkara et al. 2024). However, the star formation window spans >1 Gyr with the source only quenching within ~300 Myr of the observation time. This can be compared to >1 Gyr old stellar population observed for ZF-UDS-7329, through the direct detection of the D4000 Å spectral feature.

Open questions remain on the pathway to quiescence. Our observations have shown that, even with ground-based K-band selections assisted with >3 μm Spitzer imaging, we are able to photometrically select the oldest massive quiescent sources in the Universe (K. Glazebrook et al. 2024; T. Nanayakkara et al. 2024). One of the common themes between most of the reconstructed SFHs of z > 3 massive quiescent galaxies across various surveys is the intense SFR at early time. This is tied to the rapid buildup of stellar masses and suggests that most of these sources could have had extensive star formation in z > 5. Thus, a natural evolution of highly star-forming and obscured submillimeter galaxies at z > 5 to z ~ 3–5 massive quiescent galaxies can be expected (F. Valentino et al. 2020; A. S. Long et al. 2023). The open question is whether all z ~ 3–5 massive quiescent galaxies undergo a submillimeter galaxy phase at z > 6 or whether all submillimeter galaxies end up becoming massive quiescent galaxies by z ~ 3. Further studies linking the two populations with robust number density constraints are required to conclusively determine the evolutionary pathways of these two galaxy populations.

Color selection techniques, specially the rest-frame U − V versus V − J color distribution of galaxies have been commonly utilized to select quiescent galaxies in the extragalactic Universe (e.g., R. J. Williams et al. 2009; L. R. Spitler et al. 2014; C. M. S. Straatman et al. 2014; C. Schreiber et al. 2018a; A. C. Carnall et al. 2020). Quiescent galaxies with strong Balmer/D4000 Å breaks show redder U − V colors and bluer V − J colors. Dusty galaxies show redder V − J colors; hence, in this space, quiescent galaxies occupy a region that helps break the degeneracy of the red colors that are also observed in dusty star-forming galaxies. The quiescent quadrant can also hold dusty post-starburst galaxies, which can smooth the Balmer break of post-starburst galaxies to look like the D4000 Å break mimicking an older quiescent galaxy. JWST now opens up a new redshift and magnitude (e.g., spectroscopy of K faint sources) frontier, which eluded ground-based observations. Thus, we need to address whether our traditional definitions for quiescence still holds in this era and whether the traditional color selection techniques used in identifying quiescent systems in the z > 3 Universe are able to effectively select quiescent sources with high purity.

In Figure 18, we show the rest-frame U − V versus V − J color distribution of our galaxies. We use the C3K models fitted to the full spectrum with Prospector including emission-line contributions for this purpose. The demarcations for quiescent, red (dusty) star formation, and non-dusty (blue) star formation follows L. R. Spitler et al. (2014) criteria. Within ±0.1 mag, except for two galaxies, all of our sources fall in the quiescent region. This corresponds to an accuracy of approximately 90% for our sample and confirms that photometric selections are highly efficient in identifying quiescent galaxy candidates at z > 3.

Zoom In Zoom Out Reset image size

The dustiest source in our sample is 3D-EGS-34322, with τ_dust ~ 1.7. 3D-UDS-39102 also exhibits a higher dust content, with τ_dust ~ 1.1. Both of these sources lie in the dusty star-forming region. Our oldest quiescent galaxy, ZF-UDS-7329, shows a low level of dust, while ZF-COS-20115 has the highest dust content among the quiescent sources, with τ_dust ~ 1.5.

Galaxies with strong emission lines as identified as possible AGNs in Section 3.7 primarily fall in the quiescent region. To investigate the effects of strong lines in the UVJ classification, we show the Prospector best-fit rest-frame SEDs of our galaxies in Figure 19. ZF-UDS-8197 has the strongest emission lines in our sample, and it is clearly visible that the strong Hα and [O iii]λ5007 emission lines fall outside the V band. If these lines fell within the V band, the galaxy's position in the U − V versus V − J color space would likely shift diagonally left toward the blue, star-forming region. Since these colors are defined in the rest frame, the fact that strong optical emission lines fall outside the V band helps in more efficiently identifying quiescent galaxies photometrically, as there is less optical emission line contamination and a greater focus on the Balmer break.

Zoom In Zoom Out Reset image size

Adding effects of AGNs to Prospector fitting only has a small role for the rest-frame UVJ colors. There are four galaxies that show an absolute shift of >0.1 dex in the color–color space. There is a mean ~0.06 ± 0.06 offset in rest-frame U − V versus V − J colors for the galaxies. The maximum shift observed is 0.19 dex. None of the galaxies move away from quiescent and/or star-forming regions when AGN effects are considered in the fitting.

To investigate the effect that AGN contributions may have on rest-frame UVJ colors, we present the Prospector best-fit rest-frame SEDs used for the rest-UVJ color analysis in Figure 19. For 3D-UDS-35168, 3D-UDS-41232, and 3D-UDS-39102, the rest-UV shows a visual deviation between the two best-fit models. However, these variations are not captured by the U band. The larger deviations in spectral shape are in the rest NIR regime of the SEDs, and ~50% of our sample shows significant variations in the NIR spectral shape. For the majority of these galaxies, the shift is beyond the region covered by observed photometry. The larger observed color shifts in galaxies such as ZF-COS-20115 and 3D-EGS-18996 are driven by the increase in the rest-frame J band due to the contributions in the AGN in the SEDs. However, only the very early stages of this shift are captured by the J band. Galaxies such as 3D-EGS-34322 also show significant variation in the NIR part of the best-fit SED, when AGN effects are included in Prospector. However, this variation is not captured by the rest-frame J band; thus, the absolute shift in colors is ≲0.1 mag.

For most of our sources, the reddest F444W NIRCam band does not cover the rest-frame J band; thus, direct observational constraints are limited in constraining this turnover. We have three sources with no JWST/NIRcam imaging. For these sources, as mentioned in Section 2, we use photometric data from the ZFOURGE and 3DHST surveys. Spitzer imaging in these surveys goes beyond the rest-frame J band. Therefore, for galaxies such as ZF-COS-19589, Spitzer IRAC Channel 4 imaging provides additional constrains to the shape of the SED, albeit with relatively lower S/N, which reduces its constraining power. Based on photometric coverage, it is evident that imaging in the mid-infrared bands with instruments such as JWST/MIRI is required in z ~ 3–5 galaxies to provide tighter constraints to the contribution of AGNs to the observed SEDs of the galaxies.

To confirm whether our galaxies are truly quiescent, we use the Prospector best-fit parameters to examine the sSFRs of our sample and the distribution of the quiescent candidates within the star-forming stellar mass relation (K. G. Noeske et al. 2007). In Figure 20, we show the distribution of our sample in this space. The main locus for star-forming galaxies is defined based on the average SFR of the massive (3 × 10⁹ < M_*/M_⊙ < 10¹⁰) 3 < z < 4 galaxies in the ZFOURGE survey as defined by C. Schreiber et al. (2018a). The definition of quiescence is explored using three different methods. First, the cut at the highest sSFR is obtained with sSFR = 1/(3t_H) ~ 0.22 Gyr⁻¹, where t_H is the age of the Universe at z ~ 4. Second, the sSFR = 0.15 Gyr⁻¹ cut is utilized, which is ×10 lower than the main star-forming locus as defined above. Finally, we also use an sSFR = 0.01 Gyr⁻¹ definition, which is generally defined as the absolute cut in sSFR to obtain red and dead galaxies (e.g., G. De Lucia et al. 2024).

Zoom In Zoom Out Reset image size

All but the two galaxies that fell within the (red) star-forming region in the UVJ color selection fall within at least one of the quenching definitions. As noted by C. Schreiber et al. (2018a), this highlights that the UVJ color selection is effectively able to identify fully quenched galaxies and galaxies with very low residual star formation. We expect the massive galaxies in the latter selection are on route to be fully quenched. It is also important to highlight the distinction between our massive z ~ 3–5 quiescent galaxies and the napping galaxies identified by JWST at higher redshifts (e.g., T. J. Looser et al. 2024; V. Strait et al. 2023). Given that these sources have young stellar populations (≲100 Myr), they have considerable UV flux and therefore cannot be photometrically selected by UVJ color selection techniques. Thus, spectroscopy is crucial to discover such sources. Furthermore, these galaxies are considerably smaller (~10⁵–10⁸M_⊙), and rejuvenation at later times is expected to further build the stellar masses. Recent analysis of z > 6 star-forming galaxies shows evidence for such stochastic SFH at early times (L. Ciesla et al. 2024; R. Endsley et al. 2024; A. Pallottini & A. Ferrara 2023). Further complications arise due to spatial color variations observed within some massive galaxies (D. J. Setton et al. 2024), which hints that central regions of massive quiescent galaxies may have formed at earlier times, with residual star formation observed in the outskirts. High-resolution observations of massive and quiescent z > 3 galaxies with spatially resolved spectroscopy is required to distinguish between different quenching pathways for galaxies in the early Universe by exploring age gradients and outflow signatures (M. Park et al. 2024).

Even within the lower resolution of the Prism observations, a subsample of our galaxies shows clear NaD absorption features. In Figure 21, we show the spectra zoomed in around the region of NaD. There is no coverage for ZF-COS-20133 and 3D-EGS-27584 due to the detector gaps in NIRSpec.

Zoom In Zoom Out Reset image size

Strong NaD absorption has been observed in quiescent galaxies at z ~ 2 (S. Belli et al. 2024; M. Park et al. 2024) and in old massive quiescent galaxies at 3.2 < z < 4.6 (A. C. Carnall et al. 2024). The high EW of NaD compared to the best-fit stellar models and the blueshifted nature of the NaD feature (S. Belli et al. 2024) has been argued as evidence for AGN-driven outflows, which may have resulted in abruptly quenching the star formation. Two of our sources, ZF-UDS-7329 and ZF-UDS-6496 are observed with the NIRSpec medium-resolution gratings by A. C. Carnall et al. (2024) and show a clear enhancement of NaD compared to their best-fit model. A. C. Carnall et al. (2024) also found that ZF-UDS-7329 has an enhanced abundance of [Mg/Fe] compared to solar values with best-fit results from alf (C. Conroy & P. G. van Dokkum 2012; C. Conroy et al. 2018), suggesting an [Mg/Fe] of $0.4{2}_{-0.17}^{+0.19}$. While Na is not an α element, the dominant production mechanisms of both Na and Mg are core-collapse supernovae and asymptotic giant branch stars (C. Kobayashi et al. 2020). Thus, enhancement of Na in the ISM in galaxies can generally mirror an enhancement of α process elements in galaxies. We also note that linking α enhancements to the formation timescales of galaxies can be complicated due to the possibility of strong outflows expelling elements such as Mg, leading to a deficiency in the measured α elements (e.g., A. G. Beverage et al. 2025).

Recent cosmological hydrodynamical simulations show that AGN feedback plays a dominant role regulating star formation, leading to quenching at z > 3. Therefore, observational signatures of current and/or past effects of AGNs should be visible in observed spectroscopy. As discussed in Section 3.7, there are six galaxies in our sample with prominent broad emission lines that indicate possible AGN activity. As shown by Figure 21, a subset of our sample also shows prominent NaD absorption, including four of the six galaxies identified as possible AGNs. These two observational signatures tied together could hint at current or past AGN-driven activity in most of the galaxies in our sample. Qualitative analysis of broad AGN and NaD outflow signatures requires higher-resolution spectra than what we currently have in the Prism mode. Additionally, deeper, higher-resolution data will allow for better constraints using line ratio diagnostic analyses, such as those presented by J. A. Baldwin et al. (1981), to tighten our understanding of what powers these emission lines. Furthermore, deeper X-ray data will also add value by further constraining the nature of the AGNs in these massive quiescent galaxies at z > 3.

Ground-based spectroscopy of z ~ 3–4 massive star-forming galaxies also shows evidence for AGNs, suggesting an enhanced fraction of AGNs in higher stellar mass galaxies in the early Universe (e.g., M. Martìnez-Marìn et al. 2024). Therefore, a picture of AGN-feedback-driven rapid quenching is emerging to explain the formation mechanisms of z > 3 massive quiescent galaxies. Future mass complete spectroscopic analysis of massive z > 3 galaxies would provide tight constraints to AGN fractions and their outflow signatures in the early Universe. However, a possible complication does arise with the recent discovery of extremely compact and red objects in the early Universe (e.g., I. Labbe et al. 2025; J. E. Greene et al. 2024), which are attributed to a previously unseen population of galaxies with AGNs and possibly progenitor bulges of lower-redshift massive quiescent galaxies. The X-ray and spectral properties of these sources are different to nominal AGNs (e.g., T. T. Ananna et al. 2024; B. Wang et al. 2024); thus, deeper constraints from X-ray and rest-frame optical and NIR spectroscopy is crucial to further constrain the properties of these galaxies and investigate their role in regulating galaxy growth in the early Universe.

In the pre-JWST era, the number densities were largely constrained using photometric selection techniques. Candidates were selected with rest-UVJ color selection techniques. Deep ground-based spectroscopy was used for spectroscopic confirmations, and photometric number densities were statistically adjusted based on the contamination rate obtained by the spectroscopic observations (e.g., C. Schreiber et al. 2018a).

The parent sample of our analysis presented by C. Schreiber et al. (2018a) found a number density of 2.0 ± 0.3 × 10⁻⁵ Mpc⁻³ for 3 < z < 4 massive quiescent galaxies with a K-band magnitude cut at K < 24.5. A. C. Carnall et al. (2020) found a similarly consistent number density of 1.7 ± 0.4 × 10⁻⁵ Mpc⁻³ when similar selection criteria were applied to their photometric sample. However, when robust constrains were added to the redshift and sSFR posteriors, the number density of the A. C. Carnall et al. (2020) sample was reduced to 6.8 ± 2.4 × 10⁻⁶ Mpc⁻³.

Building upon the number densities presented by C. Schreiber et al. (2018a), we find that two further sources are classified as star-forming based on our analysis presented in Section 4.2. Both 3D-EGS-34322 and 3D-UDS-39102 did not have ground-based spectroscopic confirmations. In the initial 24 massive quiescent galaxy candidates presented by C. Schreiber et al. (2018a), five galaxies are considered to be contaminants. ZF-COS-20032 and 3D-UDS-27939 were found to be z < 2.5 dusty interlopers. Based on the MOSFIRE spectroscopic observations, C. Schreiber et al. (2018a) found that ZF-COS-20133 and ZF-UDS-8197 had rest-frame U − V versus V − J colors that moved them out of the quiescent region due to high-EW [O iii]λ5007 emission lines. 3D-EGS-18996 was found also to have moved to the star-forming region due to significant underestimation of the photometric redshift.

In Section 4.2, we recomputed rest-frame UVJ colors for the JWST/NIRspec sample. This includes all three galaxies that were ruled out to be z ~ 3–4 star-forming galaxies by C. Schreiber et al. (2018a). The rest-frame UVJ colors and the sSFRs of all three of these sources classify as quiescent. We note that we consider 3D-EGS-18996 rest UVJ quiescent due to it falling within <0.1 mag from the quiescent region in the rest-frame U − V versus V − J color space. 3D-EGS-34322 and 3D-UDS-39102 are star-forming sources and can be removed as contaminants. While ZF-COS-10559 is quiescent, it is at z = 4.3, which is beyond the z ~ 3–4 region probed by C. Schreiber et al. (2018a). Thus, once we also remove this source, we end up with the same 80% purity for quiescent sources and a number density of (1.4 ± 0.3) × 10⁻⁵ Mpc⁻³ as C. Schreiber et al. (2018a). Thus, our analysis highlights that pre-JWST number densities of massive z ~ 3–4 quiescent galaxies are largely still consistent with post-JWST results. Current cosmological and semianalytical models are able to reproduce these number densities at z ~ 3 (C. D. P. Lagos et al. 2024; F. Valentino et al. 2023) but challenges lie at z > 4 (E. J. Weller et al. 2025).