The abundance and nature of high-redshift quiescent galaxies from JADES spectroscopy and the FLAMINGO simulations

ABSTRACT

We use NIRSpec/MSA (Micro Shutter Assembly) spectroscopy and NIRCam (Near-Infrared Camera) imaging to study a sample of 18 massive (⁠|$\log M_\star /\mathrm{M}_\odot \gt 10$| dex), central quiescent galaxies at |$2\le z \le 5$| in the GOODS (Great Observatories Origins Deep Survey) fields, to investigate their number density, star formation histories, quenching time-scales, and incidence of active galactic nuclei (AGN). The data depth reaches |$\log M_\star /\mathrm{M}_\odot \approx 9$| dex, yet the least-massive central quiescent galaxy found has |$\log M_\star /\mathrm{M}_\odot \gt 10$| dex, suggesting that quenching is regulated by a physical quantity that scales with |$M_\star$|⁠. With spectroscopy, we assess the completeness and purity of photometric samples, finding number densities 10 times higher than predicted by galaxy formation models, confirming earlier photometric studies. We compare our number densities to predictions from FLAMINGO (Full-Hydro Large-scale Structure Simulations with All-sky Mapping for the Interpretation of Next Generation Observations), the largest box full-hydro-simulation suite to date. We rule-out cosmic variance at the 3|$\sigma$| level, providing spectroscopic confirmation that galaxy formation models do not match observations at |$z>3$|⁠. Using FLAMINGO, we find that the vast majority of quiescent galaxies’ stars formed in situ, with these galaxies not having undergone multiple major dry mergers. This is in agreement with the compact observed size of these systems and suggests that major mergers are not a viable channel for quenching most massive galaxies. Several of our observed galaxies are old, with four displaying 4000 Å breaks with formation and quenching redshifts of |$z\ge 8$| and |$\ge 6$|⁠. Using tracers, we find that eight galaxies host AGN, including old systems, suggesting a high AGN duty cycle with a continuing trickle of gas to fuel accretion.

1 INTRODUCTION

While precise boundaries may vary (e.g. Marchesini et al. 2023), massive quiescent galaxies are generally defined as having stellar masses |$M_\star \gt 10^{10}~\mathrm{M}_\odot $| (e.g. Glazebrook et al. 2017; Carnall et al. 2023a) and negligible star formation in the last 100 Myr (⁠|$\mathrm{sSFR_{100\, Myr}}\, {\le }\, 0.2/t_{\rm obs}$|⁠, where |$t_\mathrm{obs}$| is the age of the Universe at the redshift of the galaxy; e.g. Gallazzi et al. 2014; Carnall et al. 2024). These systems are one of the key puzzling results of the first two years of JWST operations (e.g. Valentino et al. 2023; Carnall et al. 2023a, 2024; Long et al. 2024). Their number density and physical properties are a crucial test for our models of galaxy formation (e.g. Beckmann et al. 2017; Donnari et al. 2021a, b; Hartley et al. 2023; Remus & Kimmig 2025), while the earliest and most massive systems may present a challenge even for cosmological models (e.g. Boylan-Kolchin 2023; Glazebrook et al. 2024).

On the cosmological side, the most massive and earliest quiescent galaxies are important tests of Lambda-cold dark matter (⁠|$\Lambda$|CDM) cosmology. The possible tension arises due to the early Universe lacking dark-matter haloes sufficiently massive to host these galaxies.

This appears to be the case even with the maximum possible baryon conversion efficiency of 100 per cent (Carnall et al. 2024; Glazebrook et al. 2024; Park et al. 2024; de Graaff et al. 2025, but see also Turner et al. 2025, for the effect of the chosen priors). The advantage of studying quiescent galaxies is that their stellar masses are more accurate than for star-forming galaxies, as shown by the reduced scatter of quiescent versus star-forming galaxies in scaling relations (e.g. the stellar mass plane, Hyde & Bernardi 2009; de Graaff et al. 2021). However, there are still some possible explanations for these massive quiescent galaxies that could resolve the apparent tension, including small systematic errors in the stellar mass measurements due, e.g. to non-solar abundance ratios (e.g. Beverage et al. 2021, 2025), varying initial mass functions (IMFs, e.g. van Dokkum & Conroy 2024), possible contributions from active galactic nuclei (AGN) to the continuum (Inayoshi & Maiolino 2025), and cosmic variance, making these single objects exceptions and not the rule (e.g. Valentino et al. 2023; Carnall et al. 2023b; Graaff et al. 2025; Jespersen et al. 2025).

As for galaxy formation models, massive quiescent galaxies are fascinating because some of them seem to have formed most of their stars extremely rapidly (e.g. Belli, Newman & Ellis 2019; Park et al. 2023, 2024; Beverage et al. 2025). Furthermore, these galaxies have ceased forming new stars despite their high accretion rates of cosmic gas (implied by the connection between high stellar mass and high halo mass). This disconnect between the halo accretion rate and the galaxy star formation rate (SFR) requires some physical mechanism ‘quenching’ star formation. Whatever processes have quenched massive quiescent galaxies, they appear to operate not only in the very early Universe (with galaxies quenched as early as z |$\approx 7\!-\!9$|⁠; Glazebrook et al. 2024; Weibel et al. 2024b), but also in incredibly little time (even accounting for the short time-scale set by the star-forming main sequence, SFMS, at |$z=5\!-\!7$|⁠; e.g. Clarke et al. 2024; Koprowski et al. 2024).

Many mechanisms have been proposed to quench massive central galaxies for hundreds of Myr (e.g. Man & Belli 2018; Curtis-Lake et al. 2023), but theoretical models and energy arguments seem to require some form of feedback from accreting supermassive black holes (SMBHs), through episodes of AGN (Silk & Rees 1998; Binney 2004; Di Matteo, Springel & Hernquist 2005; Bower et al. 2006; Croton et al. 2006).

This AGN feedback could act either by removing the star-forming gas through outflows (ejective feedback; Di Matteo et al. 2005; Maiolino et al. 2012; Carniani et al. 2016), or by stopping the accretion of cold gas by heating the circumgalactic medium (preventive feedback; Bower et al. 2006; Croton et al. 2006).

In support of preventive feedback, evidence from redshifts |$z=0\!-\!2$| suggests that the best predictor of quiescence is SMBH mass (Terrazas et al. 2016; Piotrowska et al. 2022) or, in the absence of direct measurements of this mass, its closest proxies, such as central velocity dispersion or bulge-to-total mass ratio (Bluck et al. 2022; Bluck, Piotrowska & Maiolino 2023). Following these findings, quenching has been attributed to time-integrated AGN feedback because (most of the) SMBH mass arises from the black-hole accretion history. This time-integrated feedback would act to prevent cosmic gas accretion onto the galaxy, which is then ‘starved’ of fuel for star formation (Peng, Maiolino & Cochrane 2015; Trussler et al. 2020; Baker et al. 2024). Most notably, these observational results match the predictions of both semi-analytical models (Bluck et al. 2022) and cosmological simulations (Piotrowska et al. 2022), which in turn correctly predict the local number density of quiescent galaxies.

However, JWST has challenged this successful picture; the discovery of a surprising abundance of massive quiescent galaxies at redshifts |$z=3\!-\!4$| (Carnall et al. 2023a), seems to defy the predictions of cosmological simulations (and semi-analytic models, SAMs, particularly at |$z=4\!-\!5$|⁠, Lagos et al. 2025), even accounting for cosmic variance in the observations (Valentino et al. 2023). One limitation of past studies, however, has been the size of the cosmological simulations used to obtain the predicted number densities. Observed densities of order |$10^{-5}$| cMpc|$^{-3}$| (Valentino et al. 2023) mean that only a few objects are found per 100-cMpc simulation box,¹ meaning that cosmic variance could still explain the observed discrepancy between observations and theory, without any change to our models.

Alternatively, assuming the mismatch between observations and theory to be representative, explaining the overabundance of observed quiescent galaxies would require a stronger or more efficient quenching mechanism at the earliest epochs, possibly some form of AGN-driven ejective feedback (e.g. Xie et al. 2024). Evidence in this last direction seems to come from the discovery of AGN-driven neutral-gas outflows in at least one-third of massive quiescent galaxies at |$z=1.5\!-\!2.5$| (Davies et al. 2024). These outflows have high mass loading, comparable to or exceeding their galaxy SFR (D’Eugenio et al. 2023; Belli et al. 2024), and their high incidence may point to a high duty cycle of AGN feedback.

One key difficulty in establishing the causes of quenching is that quenching mechanisms may have short time-scales relative to the duration of quiescence. This is particularly true for AGN feedback, because SMBHs are believed to quench galaxies over many AGN episodes, each lasting shorter than 1–10 Myr (e.g. Harrison et al. 2018; Scholtz et al. 2018, 2021; Harrison & Ramos Almeida 2024). In contrast, galaxies can stay quiescent for hundreds of Myr (e.g. Thomas et al. 2005; McDermid et al. 2015). There has been a lot of work in studying massive quiescent galaxies at |$z < 1\!-\!2$|⁠, to investigate the causes of quenching (e.g. Thomas et al. 2005; Cappellari et al. 2011; McDermid et al. 2015; Newman et al. 2018; Belli et al. 2019; Barone et al. 2022; Bluck et al. 2022; Piotrowska et al. 2022; Woodrum et al. 2024). However, by studying quiescent galaxies at increasingly higher redshifts, the younger age of the Universe necessarily leaves less room for SMBHs to be in inactive phases because the multiple AGN episodes necessary to quench the galaxy must occur more frequently in the limited time frame available. Therefore, an interesting question regarding quenching mechanisms is how many high-z quiescent galaxies are found to host AGN.

Indeed, a surprising result of the last two years has been the large number of high-z AGN candidates detected, including both Type 1 (Harikane et al. 2023; Kokorev et al. 2023; Greene et al. 2024; Matthee et al. 2024; Maiolino et al. 2024) and Type 2 (Scholtz et al. 2023; Tacchella et al. 2024), many of which appear to be overmassive compared to the host galaxy (Übler et al. 2023; Furtak et al. 2024; Mezcua et al. 2024; Juodžbalis et al. 2024b; Maiolino et al.2024) based on local scaling relations (but see Sun et al. 2025 for a different view). Their overmassive nature, if confirmed, almost certainly implies that a higher-than-expected amount of energy has already been released into the host galaxy and/or halo, and this release happened in a short amount of time. This suggests that these black holes have had significant accretion rates within their history, meaning that AGN feedback is likely to have affected the host galaxy.

However, establishing whether a galaxy is an AGN host can prove challenging (e.g. Hickox & Alexander 2018; Scholtz et al. 2023), and may require a combination of different tracers (e.g. Calabrò et al. 2023), including X-rays, emission-line diagnostics, mid-infrared, and radio emission.

A different but related problem is the lower mass boundary of massive quiescent galaxies. The discovery of low-mass, ‘mini’ or rapidly quenched galaxies at |$z=5\!-\!7$| (⁠|$M_\star \sim 10^7\!-\!10^9~\mathrm{M}_\odot $|⁠; Looser et al. 2023, 2024; Strait et al. 2023; Baker et al. 2025b) raises the question of what is the mass threshold where ‘bursty’ star formation and short-lived quenching are supplanted by long-term quenching, i.e. quiescence.

Many recent studies have focused on photometric samples (Merlin et al. 2019; Valentino et al. 2023; Carnall et al. 2023a; Alberts et al. 2024; Long et al. 2024; Ji et al. 2024a), deep spectroscopic samples at redshifts |$z=1\!-\!3$| (Park et al. 2024; Beverage et al. 2025), small samples (or even individual objects) confirmed by spectroscopy at |$z>4$| (Carnall et al. 2023b; Weibel et al. 2024a; Graaff et al. 2025; Wu 2025), and shallow spectroscopic samples at |$z>3$| (Nanayakkara et al. 2024, 2025). Each of these approaches has some disadvantages. Photometric studies may suffer from inconsistent filter coverage and are insensitive to the presence of optical AGN; deep spectroscopy at |$z=1\!-\!3$| targets the epoch when the mismatch between theory and observations is relatively small. Single-object studies may not be representative and shallow spectroscopy is not sensitive to the lower mass end of the population.

In this work, we use deep low-resolution JWST spectroscopic data to estimate an accurate number density of quiescent galaxies and to measure their lower mass envelope. We complement this sample with medium-resolution rest-frame optical spectroscopy and X-ray to radio observations to measure the fraction of AGN hosts. Finally, we combine our number density measurements with theoretical predictions from FLAMINGO (Full-Hydro Large-scale Structure Simulations with All-sky Mapping for the Interpretation of Next Generation Observations), the largest-box cosmological simulations to date, to provide a comparison of observations and theory over a meaningful survey volume.

Throughout this work, we assume the standard flat |$\Lambda$|CDM cosmology from Planck Collaboration VI (2020). Emission lines are given using vacuum (air) wavelengths in the UV (optical–NIR) rest frames (but in the analysis we use vacuum wavelengths as appropriate). We assume a Chabrier (2003) IMF.

2 DATA

2.1 NIRSpec data

In this work, we use JWST/NIRSpec multi-object spectroscopy (Jakobsen et al. 2022) from the GOODS fields (Great Observatories Origins Deep Survey; Giavalisco et al. 2004). These data were obtained as part of three programmes: JWST Programme ID (PID) GO 2198 (Barrufet et al. 2021); PID 1207 (PI: Rieke); and JADES (the JWST Advanced Deep Extragalactic Survey; PIDs 1210, 1180, 1181, 1286, 1287, and 3215; Eisenstein et al. 2023a, b). In the GOODS fields, we can combine the accurate total fluxes of deep NIRCam (Near-Infrared Camera) photometry with the superior precision of MSA NIRSpec spectroscopy (Micro Shutter Assembly; Ferruit et al. 2022). An extensive description of the NIRSpec observations and data reduction process has been provided in Bunker et al. (2024) and D’Eugenio et al. (2025). Here, we provide only a short summary, highlighting the differences compared to the publicly available data products.

All JADES NIRSpec observations use slitlets consisting of three shutters. Each observation uses a three-nod observing pattern and a variable number of dithers, as described in D’Eugenio et al. (2025). The three-nod pattern is normally used for accurate background subtraction, but can cause wavelength-dependent self-subtraction of extended sources (D’Eugenio et al. 2025). To avoid this problem, where possible, in this work, we always use the ‘two-nod’ background subtraction, which only subtracts observations that are 2 micro-shutters apart. This costs one third of the exposure time. An exception to this is the four galaxies drawn from PID 2198; these spectra were obtained from the DAWN JWST Archive² (DJA; Heintz et al. 2025), and were reduced using a combination of the jwst pipeline and msaexp.³ While these four sources do not use the two-nod background subtraction, they are all very compact (as expected from |$z>4$| quiescent galaxies; Ji et al. 2024a). The integration times range from 2 h (PIDs 1180 and 1181) to 50 h (PID 3215), with the minimum duration considered sufficient to measure reliable star formation histories (SFHs; e.g. Glazebrook et al. 2024; Nanayakkara et al. 2024).

The NIRSpec data undergo an initial processing through the JWST Calibration Pipeline applying standard detector-level corrections such as dark subtraction, linearity calibration, and cosmic ray identification. The pipeline also performs flat-fielding, background subtraction, and wavelength calibration.

For each NIRSpec exposure, we utilize the JWST grating equation models to determine the trace and wavelength solution for each source in the field of view. We then performed optimized spectral extraction, accounting for the curved traces and varying spatial profiles of the sources along the dispersion direction. Our extraction routines are optimized to handle differential slit-losses for point-like sources. To correct for these issues in the slit-loss correction and, possibly, in the overall flux calibration, we recalibrate the spectra using NIRCam photometry.

A critical component of our NIRSpec reduction process is the determination of accurate spectroscopic redshifts for the observed sources. We use an initial redshift estimate using bagpipes (Carnall et al. 2018), then refine this estimate with a combination of visual inspection, full spectrum fitting with ppxf (Cappellari 2023), and further emission-line fitting in the medium-resolution data (when available, and when emission lines are detected). The procedure is described in D’Eugenio et al. (2025).

All galaxies that we explore spectroscopically have spectroscopy obtained with the prism disperser, from PID 2198 and from JADES. These data have a spectral resolution |$R=30\!-\!300$| in the wavelength range 0.6–5.3 |$\mu$|m (Jakobsen et al. 2022). For a subset of galaxies, we also have medium-resolution spectroscopy. A fraction of these are from PID 1207, which uses a combination of G140M/F100LP and G235M/F170LP, spanning 1–1.6 |$\mu$|m with |$R=700-1200$|⁠. The remainder are from JADES which uses G140M/F070LP, G235M/F170LP, and G395M/F290LP, covering the same spectral range as the prism. Most of the medium-resolution observations are too shallow to constrain the stellar continuum (2 h). Moreover, the mask design of the JADES survey allows for spectral overlaps in the medium-grating observations (Bunker et al. 2024; D’Eugenio et al. 2025, except for the highest priority targets, e.g. Bunker et al. 2024). Although overlaps are easy to identify and remove for emission-line science (D’Eugenio et al. 2025), they may severely affect the continuum. For this reason, in this work we focus on prism spectroscopy. However, medium-resolution data are still crucial for detecting optical AGN via emission-line diagnostics. For this reason, we use these data to measure the flux of nebular emission lines, given that the spectral resolution of the gratings enables us to spectrally resolve the [N ii]|$\lambda \lambda$|6543, 6584–H|$\alpha$| complex; Section 3.2). For one galaxy from PID 1207, we also use the grating data to measure the SFH and to assess the reliability of the prism-based inference (Appendix  F). We find that the two SFHS are are consistent to within 1|$\sigma$| 1. For this comparison, we do not use the JADES gratings observations, because the JADES masks are designed by allowing spectral overlaps, which affect the continuum and can therefore bias the SFH.

2.2 NIRCam data

The NIRCam imaging data from the JWST obtained by JADES (PIDs: 1180, 1181, 1210, 1286, 3215; Rieke et al. 2023; Eisenstein et al. 2023a, b) and the JWST Extragalactic Medium -band Survey (JEMS, PID: 1963; Williams et al. 2023) are processed through the JADES custom reduction pipeline.

We start with the JWST Calibration Pipeline version 1.9.6, incorporating the latest Calibration Reference Data including dark frames, distortion maps, bad pixel masks, read noise, superbias, and flat-field files (Rieke et al. 2023).

At the detector level, we perform group-by-group corrections that include dark subtraction, linearity calibration, and cosmic ray identification, followed by ramp fitting to recover count-rate images (Rieke et al. 2023). Further calibration steps apply flat-fielding, photometric calibration, and background subtraction to ramp-fitted images (Rieke et al. 2023). We use a custom version of TweakReg to align exposures, matching sources across images to calculate relative and absolute astrometric solutions. The calibrated aligned exposures are then combined into mosaics using the JWST Stage 3 pipeline (Rieke et al. 2023).

We used the standard photometric JADES catalogues for the photometry (Rieke et al. 2023; Eisenstein et al. 2023b; D’Eugenio et al. 2025). A brief outline of these are as follows. We use the latest version of SExtractor (Bertin & Arnouts 1996) to detect sources in mosaicked NIRCam images. The detection parameters are tuned to optimize sensitivity to faint sources while minimizing spurious detections. Isophotal photometry with careful deblending allows us to measure accurate colours for blended objects. Astrometric alignment is performed using common sources across data sets to enable consistent multiwavelength photometry.

We used a combination of Kron convolved and CIRC5 photometry (which corresponds to an aperture radius of 0.35 arcsec.). We checked both and found no significant difference in the derived photometry, just differences in the goodness of fit due to extended sources and near neighbours. In each case, we chose the better fitting photometry.

2.3 FLAMINGO simulation

For comparisons with models and to help interpret our observational results, we adopt predictions from the FLAMINGO (Kugel et al. 2023; Schaye et al. 2023) project. The FLAMINGO suites are state-of-the-art cosmological hydrodynamical simulations with large volumes of up to |$(2.8\, \rm cGpc)^3$|⁠. The simulations are run using the code swift (Smoothed-particle-hydrodynamics With Inter-dependent Fine-grained Tasking, Schaller et al. 2024).

The main advantage of using FLAMINGO for this kind of study is clear: the massive objects of interest are very rare, and combined with the relatively narrow area probed by deep surveys including JWST, the cosmic variance is expected to be significant, going up to or beyond an order of magnitude (Lim et al. 2024; Xiao et al. 2024). Given the comoving volume probed by our data or other similar studies of |$10^5{-}10^6\, {\rm cMpc}^3$|⁠, only simulations larger than a cGpc size allow the probe of cosmic variance safely up to 3|$\sigma$|⁠, which is the minimum for evidencing tensions.

The suites have variations in volume and resolution as well as physics models and cosmology. For our analysis, we use its |$(1\, \rm cGpc)^3$| flagship run, L1_m8, where the mass resolution of the initial baryonic elements is |$m_{\rm gas}= 1.34\times 10^8\, {\rm \mathrm{M}_\odot }$|⁠, and the Plummer-equivalent comoving gravitational softening length is 11.2 ckpc. The reason we chose the smaller-box higher resolution simulation is that the |$1\, {\rm cGpc}$| box is already large enough, thousands of times bigger than the survey volume (see Section 5.1). This large box size enables us to explore the cosmic variance up to 3|$\sigma$| and achieve our science goals. Also, the SFH predicted from the simulations is found to be largely dependent on the numerical resolution, with the lower resolution FLAMINGO yielding a much younger SFH (and essentially zero, falling below the resolution, at |$z\, {\gtrsim }\, 7$|⁠). Although the relatively poor resolution of even the highest resolution FLAMINGO run is potentially a limitation for comparisons, increasing resolution typically leads to a lower quenched fraction (e.g. Schaye et al. 2023), because more star-forming small-scale gas is resolved, and because fewer galaxies fall below the resolution to be assigned zero SFR. Therefore, our results from FLAMINGO can be rather taken as a conservative limit for comparisons of e.g. the number densities, concerning the resolution.

The free parameters for the subgrid models employed in FLAMINGO were calibrated by machine learning to both the stellar mass function (SMF) and the cluster gas fraction from |$z\, {\lesssim }\, 0.3$| observations. The simulation employs the cosmological models from the Dark Energy Survey year three results (Abbott et al. 2022), which is a flat universe constrained by 3|$\times$|2pt correlation functions. Although this is different from the Planck cosmology assumed elsewhere in this paper, we checked, with a lower resolution FLAMINGO that adopted the Planck cosmology, that the difference is insignificant. Finally, the simulation assumes a Chabrier (2003) IMF, the same as for the observations.

3 ANALYSIS

3.1 SED fitting

We use the flexible prospector code (Johnson et al. 2021) to perform a simultaneous SED modelling of our multiwavelength photometric and spectroscopic data sets from JWST. prospector is a Bayesian Spectral Energy Distribution (SED) modelling framework designed to generate realistic galaxy SEDs from flexible stellar population synthesis models and constrain the corresponding model parameters given observational data. Source-frame galaxy SEDs are constructed using the fsps software (Conroy, Gunn & White 2009; Conroy & Gunn 2010), allowing for user-specified parametric forms or non-parametric representations of the SFH, metallicity evolution, dust attenuation, and other components.

We used the updated Medium-resolution Isaac Newton Telescope library of empirical spectra (MILES) spectral library (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011) and Mesa Isochrones and Stellar Tracks (MIST) isochrones (Choi et al. 2016). Note that these assume fixed (solar) abundances, which are not representative of the observed properties of massive quiescent galaxies at |$z\lesssim 2$| (Kriek et al. 2016; Beverage et al. 2021, 2025), and may result in overpredicted stellar ages. We assume a Chabrier (2003) IMF and a non-parametric SFH, described by a constant SFR over a set of predetermined intervals in lookback time. The model is parametrized by the logarithm of the SFR ratio between adjacent time intervals. We adopt the ‘continuity’ probability prior on this SFH model (Leja et al. 2019), which consists of a Student’s t distribution with mean 0, and standard deviation 0.3. The mean of 0 centres the prior on a constant SFR with cosmic time, while the broad wings of the t distribution allow for variations in the SFR over time, albeit with a modest penalization.

Dust attenuation assumes the Calzetti et al. (2000) attenuation law.

We perform a full spectrophotometric fitting by combining the prism spectra with multiband photometry, enabling us to incorporate information from both of these sources. We used a second-order polynomial to account for offsets between the photometry and the spectra due to wavelength-dependent slit losses. This has the effect of modelling both the spectrum and the photometry, particularly important in the case of larger extended sources where the spectrum only probes a region of the overall galaxy. We note, however, that in this process we are still fitting the photometry combined with the spectrum (so we include the shape and normalization of the photometry etc.) as this is crucial for determining the total galaxy properties. We are not simply upscaling the spectrum. This is shown in Appendix  E and Fig. E1. In the case of galaxies 12619, 171147, and 72 396, we found we required a 10th-order polynomial to better fit the observed photometry. We also test using a zero-order polynomial, i.e. not accounting for the wavelength-dependent slit losses, and find it leads to significantly poorer fits.

We test the effects of using a bursty continuity prior on the SFH (as in Tacchella et al. 2022) – this is enacted by increasing the width of the Student’s t distribution prior between the SFH bins. This has the effect of allowing greater variation in the SFR between adjacent bins, i.e. more ‘burstiness’ in the SFH. We find that this does not significantly change our results.

We also explore different priors on the stellar metallicity. This is important for two reasons. In general, for the youngest quiescent galaxies in our sample (i.e. spectral ‘post-starburst’ galaxies, e.g. Goto 2005; D’Eugenio et al. 2020; Wu et al. 2020; Suess et al. 2022), metal absorption lines have low equivalent width (EW), thus requiring high continuum signal-to-noise ratio (SNR) to reliably measure metallicity. In addition, because the spectral resolution of our prism data is insufficient to capture variations in metal-absorption features, we do not expect to constrain the stellar metallicity (Graaff et al. 2025). We start with a uniform prior in |$\log Z$|⁠, with |$\rm \log (Z/\mathrm{Z}_\odot)\in (-2,0)$|⁠. We then vary this to a higher metallicity prior to |$\rm log(Z/\mathrm{Z}_\odot)\in (-1,0.19)$|⁠, and we also test a uniform prior in linear Z. We find that whilst obtaining a multitude of different, relatively unconstrained metallicities for these galaxies, this does not alter the other stellar population properties such as formation and quenching times.

In our fiducial run, we also mask the NIRSpec spectrum for rest-frame wavelengths bluer than 3500 Å, to avoid deriving metallicity information from the relatively poorly known rest-frame UV spectrum, which may additionally suffer from strong non-stellar features, such as interstellar medium absorption (e.g. Calzetti et al. 2000) or AGN continuum emission (e.g. Carnall et al. 2023b).

The high-quality spectroscopic redshifts from NIRSpec are incorporated as priors, reducing the number of degeneracies in the SED modelling. Stellar-population properties, including ages, metallicities, and SFHs, are derived from the combined spectrophotometric constraints, leveraging the continuity SFH to obviate the lack of constraining power in the data for the earliest stages of star formation. We marginalize over the emission lines because we cannot accurately model those galaxies that contain AGN.

Fig. 1 shows an example of our prospector fits (in this case to the highest redshift quiescent galaxies found in this work, see Section 4). The observed spectrum and errors are given in grey and the observed photometry is represented by the yellow points. The best-fitting spectrum is in red and the best-fitting photometry are the black squares. The lower panel shows the |$\chi$| values for the fit of the spectra. We note that the |$\chi$| values do not take into account the marginalization of the emission lines hence the divergence seen in those regions. In some galaxies, there is an offset between the spectrum and the photometry, but the spectrum has been scaled to the photometry with a second-order polynomial (see above).

3.2 Emission-line fluxes

The emission-line fluxes for H|$\beta$|⁠, [O iii]|$\lambda$|5007, H|$\alpha$|⁠, and [N ii]|$\lambda$|6583 were measured by simultaneously fitting the continuum and emission lines with ppxf (Cappellari 2023). We used the same setup as in D’Eugenio et al. (2025), for both the prism and medium-resolution spectra.

When measuring H|$\beta$| and [O iii]|$\lambda$|5007 we always use the prism data, because prism spectra are available for all galaxies, and because they have a higher SNR than the shorter exposure grating spectra. We use the latter to measure H|$\alpha$| and [N ii]|$\lambda$|6583, given that the medium-resolution data has sufficient spectral resolution to separate H|$\alpha$| from the [N ii]|$\lambda \lambda$|6548,6583 doublet. Measurement uncertainties are obtained by adding random noise to the observed spectra and repeating the fit. Random noise is derived by using the pipeline uncertainties. We tested sampling from the best-fitting residuals in a neighbourhood of each spectral pixel (e.g. Looser et al. 2023), and found similar results. This procedure also varies the shape of the underlying stellar continuum, but artificially increases noise (by adding random noise to already noisy data). We reduce the uncertainties derived by a factor |$\sqrt{2}$| to take into account the additional noise introduced by this procedure.

3.3 Morphological fitting

To fit the morphology of the quiescent galaxy sample, we use the code PySersic (Pasha & Miller 2023). This software forward models and fits the light distribution of the galaxy assuming a parametric model. We use a standard, axisymmetric (Sersic 1968) profile. For the point spread function (PSF), we use model PSFs from the JADES collaboration as in Ji et al. (2024a, c). As input, we use the F444W images because this is the only wide-band filter available for all targets, whilst also not being contaminated by strong emission lines from AGN (we also test using F200W to make better use of the SW PSF and find it makes no significant difference to our results). For our quiescent galaxy sample, this corresponds to rest-frame optical continuum emission from the stellar populations, given that at |$z < 5$| we can rule out strong emission line in this filter and a major contribution from hot-dust emission due to AGN. We create 3 arcsec|$\times$| 3 arcsec square cutouts centred on each target galaxy; this size is large enough to provide source-free regions to estimate any residual background, but small enough to save computational time. We use Source Extractor to model any remaining background in the cutouts and to perform source detection. We fit a single Sérsic (1968) profile to every source in the cutout, except those with fluxes below 10 nJy, which we model as point-source components due to their small size and low SNR. In almost all cases, the primary source contains almost the entire flux in the image, but we still fit the remaining sources to avoid biasing the model by trying to account for excess light away from the central source. We use NUTS (No U Turn sampling) to sample the posterior in order to determine the best-fitting model and to obtain more realistic uncertainties. The best-fitting half-light radii and Sérsic indices are reported in Table B1.

4 SAMPLE

4.1 JADES sample

We selected our samples from spectroscopic data in the GOODS-S and GOODS-N fields. We limit the selection to galaxies with prism data, and to two programmes: JADES and PID: 2198 (see Section 2.1). As we are interested in higher redshift quiescent galaxies we imposed a redshift cut of z |$\ge$| 2.

First, we select candidate quiescent galaxies based on UVJ colours obtained by template fitting photometry with the fitting code eazy (Brammer, van Dokkum & Coppi 2008), as described in Hainline et al. (2024). eazy uses galaxy templates to perform a grid search based upon redshift. The templates used are a combination of the eazy ‘v1.3’ templates (Brammer et al. 2008; Erb et al. 2010), alongside mock templates from the Jades Extragalactic Ultradeep Artificial Realizations (JAGUAR) simulations (Williams et al. 2018). These are supplemented by new templates designed using fsps (Conroy et al. 2009) which were added to those of Coe et al. (2006).⁴

We use the UVJ criteria of Schreiber et al. (2015) with the fast quenching addition of Belli et al. (2019). Specifically, the dividing line given by

where U, V, and J are rest-frame magnitudes in the same filters as Williams et al. (2009). However, as we are looking for rare, high-redshift quiescent galaxies, we want to ensure that we recover all quiescent galaxies with spectra. Therefore, we include the 1|$\sigma$| uncertainties, thereby selecting any candidate galaxy that falls into the extended Schreiber et al. (2015) criteria (Belli et al. 2019) if the U, V, J colours are perturbed by 1|$\sigma$|⁠. Also, as we are looking for massive quiescent galaxies, we impose a stellar mass cut of |$M_\star \gt 10^{10}\,\,\mathrm{M}_\odot $| based on the best-fitting eazy masses. We test our best-fitting prospector masses vs the template derived eazy masses finding no systematic offset and a scatter of within 0.2 dex.

Next, to select the final sample from these candidates, we imposed the time-dependent criterion based on the specific SFR (sSFR) obtained via SED fitting (Section 3.1).⁵ This criterion is given by

where |$\rm sSFR_{100Myr}$| is the SFR averaged over 100 Myr divided by the total formed stellar mass by the time of observation, and |$t_\mathrm{obs}$| is the age of the Universe at the redshift of the galaxy. This criterion has been used previously in many studies and is generally considered equivalent to a colour cut in the UVJ colour–colour diagram (e.g. Gallazzi et al. 2014; Carnall et al. 2024). Interestingly, in our case we find that quiescent galaxies that fall outside of the standard UVJ quiescent region are still clearly quiescent based on the sSFR cut. We also, as a double check, visually inspect the spectra and see clear evidence for strong Balmer breaks and a lack of emission lines. This motivates further study to explore how many high-z quiescent galaxies are lost by relying on colour selection alone at high redshift where significant extrapolation can be required (e.g. Antwi-Danso et al. 2023; Gould et al. 2023; Valentino et al. 2023).

The sSFR selection results in a sample of 21 quiescent galaxies, 17 from JADES and 4 from PID2198 (which includes CANDELS GS-9209; Carnall et al. 2023b). After visual inspection, we reject three galaxies from this sample; one due to a spectrum contaminated by short circuits in NIRSpec (Rawle et al. 2022), and two that exhibit a combination of a Balmer break and dust-obscured star formation – our treatment of emission lines in the spectrum are unable to disentangle this kind of system.

This gives us a final sample of 18 massive quiescent galaxies.

The UVJ colours of the final sample are shown in Fig. 2. Schreiber et al. (2015) criteria is given by the dashed black line and the fast quenched extension (Belli et al. 2019) is given by the grey dashed line.

We find that in order to include all our quiescent galaxies based on their 50th percentile eazy colours, we would require a selection criteria of

shown as the orange-dotted line in Fig. 2. This motivates the exploration of the systematics involved in the colour selection at higher redshift.

We note that during the selection process we also identified two low-mass quiescent satellite galaxies below the mass cut. Moreover, we identified another galaxy above the mass cut that is clearly a satellite. Because we are interested only in internal quenching, we reject these three galaxies, which are clearly satellites, based on their proximity to higher mass galaxies at the same redshift (one of which has been already reported by Sandles et al. 2023). Nevertheless, the fact that our methods correctly identify even the two low-mass satellite galaxies as quiescent means that the deepest spectroscopic data we use are sensitive enough to probe the lowest mass envelope of the quiescent population. We explore these galaxies’ spectra and SFHs in more detail in Appendix  A.

For our final sample of 18 massive quiescent galaxies, we report galaxy NIRSpec IDs, coordinates, Balmer-break strengths, and D|$_n$|4000 index as well as other derived properties in Appendix Table B1.

Fig. 3 shows red-green-blue (RGB) images of all 18 massive quiescent galaxies in our sample. All appear classically quiescent (i.e. red, with no blue clumps or other signs of recent star formation) on the basis of the NIRCam imaging. We model all the quiescent galaxies and other sources in the imaging with a single-component Sérsic profile with a freely varying Sérsic index (see Sections 3.3 and 5.9).

Fig. 4 shows stellar mass versus spectroscopic redshift (left-hand panel) and SFR (averaged over 100 Myr) versuss stellar mass (right-hand panel). From the left-hand panel, we see that we have massive quiescent galaxies from redshift 2 (our lower cut) to redshift 4.7.

We find five galaxies above redshift four – two of which (1397 and 80660) have not previously been published, one of which was previously identified as a little red-dot candidate (Kokorev et al. 2024), one of which (8290) had its spectra previously published in Barrufet et al. (2025) but not explored (i.e. fit), and one (8777) whose R100 spectra were published in Barrufet et al. (2025) and whose R1000 spectra was explored in detail in Carnall et al. (2023b). This galaxy (ID 1297, CANDELS ID GS-9209) is the highest redshift quiescent galaxy in our sample at z = 4.7 and hosts a broad-line region (BLR) AGN (Carnall et al. 2023b). As mentioned above, 1397 has not been published previously, probably because of the detector gap affecting part of the prism spectrum. However, we still see clear evidence of a Balmer break, which, complemented by JADES photometry, suggests a lack of high-EW emission lines (estimating EWs with the photometry and model spectrum gives EW(H|$\beta$| + [O iii]) = |$7\pm 0.3$| Å and EW(H|$\alpha$| + [N ii]) = |$25\pm 1.3$|⁠). The two quiescent galaxies from JADES both reside in GOODS-N, at redshifts 4.1 and 4.4, respectively.

The observed spectra and photometry for these five (z > 4) galaxies can be seen in Fig. 1 along with the best-fitting models. We show the remaining spectra in Fig. 10 and in Appendix Figs C1 and C2. We see clear evidence for strong Balmer breaks and absorption features, confirming these galaxies’ quiescent nature.

We find that our galaxies lie within the stellar mass range of |$10^{10}-10^{11.3}\,\mathrm{M}_\odot $|⁠, with the majority having a lower mass than other comparable samples at this redshift (e.g. Carnall et al. 2024; Nanayakkara et al. 2024).

We next explore the galaxies’ position relative to the SFMS (Brinchmann et al. 2004; Noeske et al. 2007; Speagle et al. 2014; Baker et al. 2023; Popesso et al. 2023). This relation describes the strong correlation between the SFR of a galaxy and its stellar mass, which stems from a combination of the molecular gas main sequence, which describes the amount of gas within a galaxy relative to its stellar mass, and the Schmidt–Kennicutt relation, which encodes information on the SFR of a galaxy relative to its amount of gas (Baker et al. 2022, 2023). Although the SFMS does not directly encode a causal relation (e.g. Baker et al. 2023), it is still a useful benchmark for studying the evolutionary stage of galaxies, even when the gas content is known (Tacconi, Genzel & Sternberg 2020; Scholtz et al. 2024).

Fig. 4 (right-hand panel) shows the galaxies’ positions relative to the SFMS at redshifts 2 and 4 (grey-dashed and dotted lines from Popesso et al. 2023). We see that all our quiescent galaxies reside far below the SFMS. This is to be expected, as they are all defined as quiescent on the basis of several different criteria. This shows that over a timescale of the past 100 Myr there has been no considerable star formation activity.

4.2 FLAMINGO sample

Based on the characteristics of the observational data, we select galaxies from the FLAMINGO simulation. Specifically, we select massive quenched galaxies from FLAMINGO with |${M_\ast \, {\ge }\, 10^{10}\, {\rm M}_\odot }$| and |${\rm sSFR_{100Myr}}\, {\le }\, 0.2\, /\, t_{\rm obs}$|⁠, mimicking the observational selection. Although the stellar mass of the JADES samples extends down to |${\simeq }\, 10^{10}\, {\rm M}_\odot$|⁠, we apply this conservative cut because the completeness of the samples is not clearly known. We tested with several other selections, such as the mass cut of |${M_\ast \, {\ge }\, 10^{10.4}\, {\rm M}_\odot }$|⁠, and found that the results almost always vary within a factor of two. The sSFR is measured by summing up the mass of star particles newly born during the last 100 Myr, and by dividing it by the total formed stellar mass as well as 100 Myr. The selection of simulated galaxies is made independently of the two snapshots at |$z\, {=}\, 2.5$| and 3.5. This is because the observational data cover a wide range of redshifts, during which galaxies undergo substantial evolution in both physical properties and number density. We also add a lognormal scatter of 0.3 dex to the stellar mass and SFR originally from the simulation before the selection, to represent the observational errors in the derived quantities.

5 RESULTS

5.1 Number densities of quiescent galaxies

One of the key questions we want to address is whether the high reported number density of massive quiescent galaxies could still be consistent with models, given a sufficiently large box, and given the higher accuracy of spectroscopic SFH measurements. This enables us to better understand the population of these galaxies – providing a stronger comparison to theoretical modelling and galaxy formation and evolutionary physics. In order to properly measure the comoving number density, we need an accurate estimate of the survey volume.

However, due to the complicated selection function of JADES spectroscopy (e.g. Bunker et al. 2024; D’Eugenio et al. 2025), the actual volume surveyed by this spectroscopic study is difficult to calculate directly.

To avoid this issue, we used a photometrically selected sample (Ji et al. 2024a) obtained from products based on the NIRCam UVJ colours of galaxies in the GOODS-S JADES + FRESCO footprint (Ji, private communication). These were taken from a survey area of 77.1 arcmin|$^2$| (Eisenstein et al. 2023a, b).

We will use the volume probed by this photometric sample to compare with the simulations, after applying completeness and purity corrections derived from spectroscopy. To compute the volume, we simply take the survey area (⁠|$\Omega$| = 77.1 arcmin|$^2$|⁠), then the differences between the volume of two cones corresponding to the redshifts probed. This corresponds to the equation

where |$d_{z_2}$| is the comoving distance at the upper redshift range and |$d_{z_1}$| is the comoving distance at the lower redshift range. We divided the volume into two redshift ranges, z = 2–3, which corresponds to |$2.73\times 10^5\, {\rm cMpc}^3$| and z = 3–4.5, which corresponds to |$3.78\times 10^5\, {\rm cMpc}^3$|⁠. One of our galaxies, ID 8777, the highest redshift one, is beyond z = 4.5, thus outside the volume. The volume completeness of the samples, however, is not clearly known while it is likely to drop at the edge. We therefore just adopt the volume as described for the calculation of the number density, but found that the conclusion does not change when expanding the volume to the redshift of 8777 or even to z = 5.

To compute the photometric number densities we divide the number of quiescent galaxies found within the volume by the volume itself.

The key advantage we have compared to other photometric number density approaches is that we have spectra, so we can determine false positives and false negatives for these photometric sample. Just by exploring the UVJ colours, we expect to find galaxies that may be missed by UVJ colour selection alone (see Fig. 2).

We verified that we recovered this photometric sample with our spectroscopic approach. To do this, we match our spectroscopic sample to the photometric sample. We find that there are nine targets in common within the given areas. Eight of these are selected in our spectroscopic sample (6620, 8290, 8777, 11252, 12619, 171147, 199773, and 200733), whilst one (6887) we manually removed as the spectrum was contaminated by shorts. This means that our spectroscopic sample finds 100 per cent of the photometric sample for which we have spectra.

We next look through our spectroscopic sample to find any galaxies missing in the photometric sample. This is key as we are then able to correct our number densities for quiescent galaxies missed by photometric selection alone. We find three spectroscopically confirmed quiescent galaxies that are missing from the photometric sample (12214, 1397, and 197911). This means we apply a correction factor to our photometric number densities based on the number of quiescent galaxies we expect the photometric sample to have missed. We multiply our number densities by the fraction of the actual number of quiescent galaxies (photometric plus those missed) to the photometrically selected number. For the |$z = 2\!-\!3$| bin, this correction factor is 5/4. For the |$z= 3\!-\!4.5$| bin, this correction factor is 6/4. For the combined bin, this corresponds to 11/8.

Our observed number densities are reported in Table 1. We also report an observed number density of 1.28|$\rm \times 10^{-5}\,\,[cMpc^{-3}]$| for the z = 4–5 range, but caution that this is based on photometrically identifying two galaxies in this region and spectroscopically finding three.

We compute the 1|$\sigma$| uncertainty range for the observational number densities by dividing the survey area into 25 subareas, and thereby drawing 10 000 bootstrap samples. For each bootstrap sample, the age of each galaxy is also altered within the derived uncertainty.

To compute the model predictions of the number densities, we divide the simulation box of |$(1\, \rm cGpc)^3$| into several sub-boxes each of which has the same volume as the observation. This allows us to take into account cosmic variance, i.e. the random chance to find a volume with the observed number densities, thus helping to quantify the significance of the discrepancy with the model. The procedure results in 3375 and 2197 sub-boxes for the |$z\, {=}\, 2{-}3$| and |$3{-}4.5$| volumes, respectively. This is a large enough number of sub-boxes to explore the cosmic variance up to about |$3.5\, \sigma$|⁠.

Fig. 5 shows the number density of massive quiescent galaxies in these observed volumes compared to the number found in the FLAMINGO simulation. The cumulative number density is presented as a function of |$t_{50}$|⁠, which is defined as the lookback time by when a galaxy formed 50 per cent of its total formed stellar mass. This means that |$t_{LB}$| = 0 corresponds to our observed number densities, whilst increasing |$t_{LB}$|⁠, decreases the number of quiescent galaxies observed with that formation time. In both redshift ranges, the simulation predictions disagree with the observed number densities at about or greater than 3|$\sigma$|⁠. At higher redshifts |$z=3\!-\!4.5$|⁠, our number densities are an order of magnitude higher than the average predicted by FLAMINGO. Although comparisons with FLAMINGO have not been investigated before, a discrepancy with predictions from theory has already been reported by other photometric studies (e.g. Valentino et al. 2023; Carnall et al. 2023a), and is also seen in other numerical simulations (see following Section 5.2 and Fig. 6). However, using the large volume of FLAMINGO, it is clear that cosmic variance, by itself, is unlikely to explain the observations, with the number density at |$z=3\!-\!4.5$| close to a 4|$\sigma$| fluctuation. The use of spectroscopy rules out significant contamination of our sample as a possible explanation. The use of prism spectroscopy enables us to measure much more precise stellar ages compared to photometry alone (e.g. Nanayakkara et al. 2024).

5.2 Comparison of number densities

It is worth considering how our computed comoving number densities from the observations and the FLAMINGO simulations compare to those in the literature. Extensive comparison to large-scale cosmological simulations (e.g. eagle, Illustris-TNG, simba, Crain et al. 2015; Schaye et al. 2015; Pillepich et al. 2018; Springel et al. 2018; Davé et al. 2019) and SAMs (e.g. shark, gaea, galform, Lacey et al. 2016; Lagos et al. 2018, 2024; De Lucia et al. 2024), have already been analysed in-depth in Lagos et al. (2025). However, in that case they were limited on the observational side by the lack of spectroscopic correction (the single point from Nanayakkara et al. 2025, not withstanding) to their photometric samples so were unable to adjust for missed quiescent galaxies or interlopers. They were also limited on the simulations side by the box-size, therefore, being unable to probe or rule out cosmic variance within the simulations themselves. With the addition of our 18 strong spectroscopic sample being able to correct our photometric sample, alongside the large-scale 1 cGpc|$^3$| FLAMINGO box, we have the ability to probe these two effects.

Fig. 6 shows comoving number densities versus redshifts for massive quiescent galaxies based on observational surveys, cosmological simulations and SAMs. The red points show the number densities computed in this work with our spectroscopically corrected photometric sample. This includes the z = 2–3, 3–4.5, and 4–5 bins (see Section 5.1 for more details on this computation.) The solid blue line corresponds to the number densities from FLAMINGO with the shaded regions corresponding to 1|$\sigma$|⁠, 2|$\sigma$|⁠, and 3|$\sigma$| errors for a volume of |$3\times 10^5\, {\rm cMpc}^3$| (an average volume for our data), again testing cosmic variance. This is similar to Fig. 5, but the focus here is on the redshift evolution not the stellar ages.

We can again straight away see the discrepancy between our observed number densities and the simulated number densities. There is a 2|$\sigma$| difference for the z = 2–3 bin, but greater than 3|$\sigma$| for the higher redshift bins. This is highlighting that we likely need both earlier and stronger feedback within the simulations. We also see that for the highest redshift bins, cosmic variance within the simulations cannot account for this discrepancy.

We can also compare these number densities to results from the literature. The other points in Fig. 6 correspond to other photometric surveys including Carnall et al. (2023a), Gould et al. (2023), Valentino et al. (2023), and the single point based on spectroscopy from Nanayakkara et al. (2024). The dotted grey line corresponds to the Illustris-TNG300 cosmological simulation (Pillepich et al. 2018; Springel et al. 2018) and the dashed grey line corresponds to the prediction from the SAM shark (Lagos et al. 2018, 2024), both of which are obtained from Lagos et al. (2025).

It is immediately apparent from the figure that we see lots of variations between the different measures. Some of these are easily explainable due to simple selection effects. In our work we are selecting massive quiescent galaxies as those with |$M_\star \gt 10^{10}\,\mathrm{ M}_\odot$|⁠. Gould et al. (2023) for example select massive quiescent galaxies with |$M_\star >10^{10.6}\,\mathrm{ M}_\odot$| and in Valentino et al. (2023) they select two samples, one with |$M_\star >10^{10.6}\,\mathrm{ M}_\odot$| and the other lower mass bin with |$9.5 < M_\star < 10^{10.6}\,\mathrm{ M}_\odot$|⁠. In Fig. 6, we plot the |$M_\star \gt 10^{10.6}\,\mathrm{ M}_\odot$| reported values as this most closely replicates our sample (our lowest mass spectroscopically confirmed massive central galaxy in GOODSS has |$M_\star =10^{10.4}\,\mathrm{ M}_\odot$|⁠). Despite this, it is natural to assume the numbers in Gould et al. (2023) and Valentino et al. (2023) would be increased by switching to a |$M_\star \gt 10^{10}\,\mathrm{ M}_\odot$| selection. If we combine the number densities for both of the bins from Valentino et al. (2023) it becomes more consistent with our estimate from this work. Even with just the |$M_\star \gt 10^{10.6}\,\mathrm{ M}_\odot$| bin we are within around 1|$\sigma$|⁠. This is important as Valentino et al. (2023) have a larger area (145.1 arcmin|$^2$| compared to our 77.1 arcmin|$^2$|⁠) due to combining photometry from multiple different surveys. This is an important check that we are not being too strongly affected by cosmic variance within our observations and remain consistent with number densities calculated in larger volumes.

There appear to be larger tensions with the number densities computed in Carnall et al. (2023a), particularly in the highest redshift bin of z = 4–5. However, as was noted in Valentino et al. (2023), the The Cosmic Evolution Early Release Science Survey (CEERS) field appears to have an overdensity of massive quiescent galaxies compared to other comparable fields, showing the role cosmic variance can play and this is likely to account for the excess seen in Carnall et al. (2023a).

For the simulations and SAMs, of upmost importance are the various feedback prescriptions for galaxy quenching (e.g. Vogelsberger et al. 2020; Lagos et al. 2025). The exact implementation and methodology for this feedback directly relates to the number of massive quiescent galaxies found at high-z. In other words, the observed number densities of massive high-z quiescent galaxies are a key test of the strength and time-scales of feedback prescriptions within cosmological simulations and SAMs. Lagos et al. (2025) found that just among simulations and SAMs there are significant differences for the predicted number densities of high-z relating to these feedback implementations, with some predicting better at certain mass scales and redshift ranges.

In Fig. 6, we overplot number densities computed in Lagos et al. (2025) for the Illustris-TNG300 (Pillepich et al. 2018; Springel et al. 2018) cosmological simulations and the shark (Lagos et al. 2018, 2024) SAM. We use the number densities corresponding to quiescent galaxies with |$M_\star \gt 10^{10}\,\mathrm{ M}_\odot$|⁠. We also use the most generous sSFR threshold of sSFR|${\lesssim }\rm 10^{-10}yr^{-1}$|⁠, although note that this makes little difference.

One of the first aspects that becomes apparent is that Illustris-TNG300 dramatically underestimates the number of quiescent galaxies above z |$\sim$| 3.5, with zero quiescent galaxies above z = 4. It was noted in Lagos et al. (2025) that the actual mass threshold for selecting massive quiescent galaxies in Illustris-TNG300 was not mass dependent as all were very massive. Whereas shark is better able to match the high-z massive galaxy number densities, but (as in Lagos et al. 2025) we caution that this is because of the number of massive galaxies with |$10^{10}\,\mathrm{M}_\odot < M_\star < 10^{10.5}\,\mathrm{M}_\odot $|⁠, whereas all our spectroscopically observed galaxies with z > 4 have |$M_\star \gtrsim 10^{10.5}\,\mathrm{M}_\odot $|⁠. If we take the shark number densities for |$M_\star \gtrsim 10^{10.6}\,\mathrm{M}_\odot $| we are likely to see significant offsets to the observed number densities in this work and other studies.

FLAMINGO under predicts at all redshifts with the largest deviation at z > 3, but it does not decline beyond z = 3.5 as rapidly as Illustris-TNG300.

Overall, this motivates the needs for larger observational surveys with spectroscopically confirmed high-z quiescent galaxies, in addition to cosmological simulations with both stronger and faster acting feedback prescriptions.

5.3 SFHs

Fig. 7 shows a selection of the individual SFHs of this sample of massive quiescent galaxies as inferred from our spectra. The remaining SFHs are shown in Appendix D, Figs D1 and D2. Technically, these are a composite of the possible SFHs inferred from prospector, where the red line corresponds to the 50th percentile of the SFR distribution in each time bin, while the shaded region refers to the 16th and 84th percentiles.⁶ The formation and quenching times (⁠|$t_{50}$| and |$t_{90}$|⁠) are indicated by the dashed blue and red lines, respectively. These times are defined by when the galaxy formed 50 and 90 per cent of its total formed stellar mass (i.e. not accounting for mass loss due to stellar evolution). None of the SFHs show signs of recent activity, but this is imposed by the sSFR criteria. Some galaxies appear to have formed their stars in a single burst (e.g. 72 396), whilst others appear to have a more extended SFH (e.g. 1397). The galaxies with the strong 4000-|$\rm{\mathring{\rm A}}$| breaks show the oldest SFHs such as 11252, 12214, and 72396.

More broadly we can see that our SFHs split into two main categories – those that are more extended, with residual star formation over a longer period of time, and those that appear to be consistent with a single burst. This can be clearly seen in the differences between the formation time (dashed blue line) and the quenching time (dashed red line). Galaxies with more extended SFHs have a much greater difference between formation and quenching times, showing that they continued to build up their mass over a long epoch. Whereas for galaxies with a single burst of star formation, they have very similar formation and quenching times; all star formation activity appears to have occurred within that short period of time (i.e. |$\le$|100 Myr), with the galaxy forming hundreds of solar masses of stars per year.

None of the galaxies appear to show statistically significant evidence of rejuvenation events, although we remove by selection any substantial rejuvenation event in the last 100 Myr. Galaxies undergoing recent rejuvenation can easily wash out the Balmer break and therefore appear non-quiescent (Witten et al. 2025). In spite of this, we do see some galaxies with signs of double peaks in their SFHs, possibly suggesting short periods of reduced activity before the star formation reignites again. However, this is not statistically significant.

5.4 Formation and quenching times

Fig. 8 shows the formation redshift (⁠|$z(t_{50})$|⁠) and the quenching redshift (⁠|$z(t_{90})$|⁠) versus the observed redshift. The dashed line corresponds to the one-to-one redshift line. We see from the figure that we have a range of formation and quenching redshifts, with a significant fraction of the quiescent galaxies forming earlier than redshift 6. This means that these galaxies formed during some of the earliest stages of the Universe’s history, i.e. within the first billion years. We also find several galaxies with quenching times above redshifts of 6 suggesting that they have been quiescent for a substantial period of time. This also indicates that they must have built up their mass early on, as they are observed as massive quiescent galaxies. This is consistent with their spectral features showing evidence of 4000-|$\rm{\mathring{\rm A}}$| breaks (Section 5.5). This is consistent with previous findings of early formation and quenching times for individual or small samples of galaxies at high-z (e.g. Carnall et al. 2024; Glazebrook et al. 2024; Graaff et al. 2025) as well as the high-z quiescent galaxy found in Weibel et al. (2024a).

However, since our samples are spread across a wide range of |$z_{\rm obs}$|⁠, galaxies with similar quenching redshifts could have been quiescent for various amount of times. Therefore, it can be more appropriate to directly investigate |$t_{90}$| or |$t_{50}$|⁠, rather than the corresponding redshifts, to discuss their ages. This is shown in Fig. 9. It nicely demonstrates that a majority of those with earlier quenching redshift have been quiescent for a longer period, thus older on average. It also enables us to better visualize the limits on formation and quenching times set by the age of the Universe (as indicated by the black dashed line).

5.5 High-z 4000-Å breaks

In this section, we explore the strong 4000-Å break observed in four galaxies in our sample, as shown in Fig. 10, and why this fits with our early formation time-scales for our quiescent galaxies.

The 4000-Å break is a significant spectral feature observed in galaxies with stellar populations older than |$\approx 1\!-\!2$| Gyr. This feature, also sometimes referred to as the D4000 break from the name of an empirical spectral index (Bruzual 1983), is characterized by a discontinuity in the spectrum around 4000 Å, where there is a marked decrease in flux density on the blue side compared to the red side. However, due to the nature of the 4000-Å break, it is not a sharp discontinuity. Another key signifier of a 4000-Å break is the smaller flux drop that is typically seen around 4300 Å.

The 4000-Å break originates primarily from the cumulative effect of numerous absorption lines from neutral hydrogen and ionized metals in the atmospheres of cooler stars. The break is quantitatively defined by the ratio of the average flux density in two narrow wavelength bands on either side of 4000 Å. Specifically, it is often measured as the ratio between the average flux density in the range 3850–3950 Å and the range 4000–4100 Å (D|$_n$|4000 index; Balogh et al. 1999).

The strength of the 4000-Å break (i.e. D|$_n$|4000) serves as a robust diagnostic tool for determining the age and metallicity of stellar populations in galaxies. Older stellar populations exhibit a stronger D|$_n$|4000 break due to the higher prevalence of cool stars and increased metal blanketing. Conversely, younger stellar populations, which contain a higher fraction of hot blue stars, exhibit a weaker D|$_n$|4000 break. A D|$_n$|4000 above a value of 1.4 is also related to the metallicity of stellar populations, as well as to their age (Kauffmann et al. 2003a).

We report D|$_n$|4000 indices for the quiescent galaxy sample in Appendix Table B1. We find a range of values from 1.23 to 1.62.

Fig. 10 shows the spectra for the four galaxies in our sample that show features consistent with a strong 4000-Å break. Overplotted in blue is the galaxy from Glazebrook et al. (2024) which has been shown to have a strong 4000-Å break /bf (⁠|$D_n 4000=1.56\pm 0.18$|⁠) consistent with an old stellar population and high formation and quenching time. We can clearly see that the breaks match between the four galaxies within our sample and that of the Glazebrook et al. (2024) galaxy. The |$D_n 4000$| values measured for our four galaxies are consistent with those of Glazebrook et al. (2024). This further supports the early formation and quenching times shown in Figs 8 and 9.

5.6 In-situ versus ex-situ star formation in FLAMINGO

A key limitation of a reconstructed SFH is the inability to trace where stars formed, i.e. whether they were born from gas already in the galaxy (in situ), or if they were accreted during a merger event (ex situ). However, numerical simulations can offer insight into the likely formation pathways. In this section, we explore the in-situ versus ex-situ star formation for massive quiescent galaxies in FLAMINGO. Fig. 11 shows the mass build-up for the massive (⁠|$M_\star \, {\ge }\, 10^{10}\, {\rm M}_\odot$|⁠) and very massive (⁠|$M_\star \, {\ge }\, 10^{11}\, {\rm M}_\odot$|⁠) galaxies from the simulation at redshifts of 2.5 and 3.5. The black curve shows the total stellar mass as a function of redshift, whilst the red curve shows the amount built up in situ. Both the medians and individual curves for 50 randomly selected samples that satisfy the same quenching criteria as in Section 4.2 are shown.

The in-situ mass from the simulation is calculated as follows. First, for each simulated galaxy, we identify its main progenitors across cosmic time using the simulation merger tree. Next, for each star particle associated with the main progenitor at each given snapshot, we identify the redshift of its birth, which is available from the simulation product. Then we go to the first snapshot after the particle’s birth and see if at that time the particle belongs to the main progenitor or not. If it belongs to the main progenitor, then we identify the star particle as in situ. Due to the finite-time intervals of |${\simeq }\, 100\,$|Myr between the simulation snapshots, we cannot rule out the cases that the star particle was born ex situ but merged into the galaxy between the birth and the first snapshot after the birth. In principle, this may therefore lead to an overestimation of the in situ fraction. Given the low ex situ or merger fraction overall as indicated by the results, it is unlikely to affect the results significantly.

The key takeaway from Fig. 11 is that, on average, the in-situ mass build-up explains the vast majority of the total mass formed. In contrast, the ex-situ mass build-up via mergers appears to be contributing always less than 40 per cent to the galaxies’ growth, even though the ex-situ fraction slowly increases with galaxy mass and time. This can also be seen by the smoothness of the black median curve. In the case of frequent and dominant major dry mergers, we would expect to see an almost stepwise increase in the mass by factors of up to two, due to the almost instantaneous increase of the mass of the resulting galaxy. We do not see this from the simulation results. The simulations predict, instead, that only a minority of the galaxies in our observed sample will have undergone a major gas-poor merger event over each |${\simeq }\,$|100 Myr. In fact, the mock galaxies have undergone about one and a half (one) major dry mergers in total by |$z\, {=}\, 2.5$| (3.5) on average, where major dry mergers are defined as mergers with the stellar mass ratio greater than 10:1. This corresponds to 0.6 such mergers per Gyr. In addition, not even all these mergers contribute significantly to the stellar mass budget at |$z_{\rm obs}$|⁠, because the impact from the mergers with a long lookback time will be dominated quickly by the in-situ growth since then. Given the low rate of 0.6 Gyr⁻¹, only a negligible fraction of galaxies at a given |$z_{\rm obs}$| are expected to have had major dry mergers recently enough to possess a significant ex-situ mass budget. This is why the black and red median curves are so close to each other.

On the other hand, they are typically found to undergo minor mergers, which are much more frequent and thus contribute more significantly to the stellar mass. However, the resolution of FLAMINGO does not allow us to probe the impact of mergers with a mass ratio lower than 10:1 at |$z\, {\gtrsim }\, 7$| for objects like our samples. In addition, assessing the impact of minor mergers depends on the definition of such mergers, as well as science goals.

The case of wet major mergers is slightly more complicated. This would result in a stepwise increase in the total baryon budget, but the star formation resulting from any gas inflow would then be counted as in-situ star formation. This could very likely ‘wash-out’ any stepwise increase seen in the figure, assuming that the newly acquired gas would be consumed gradually over a long period, not instantaneously. Therefore, we caution that we cannot rule out the impact of wet major mergers in the same way as we can rule out dry major mergers (although at early epochs it is difficult to distinguish between wet mergers and gas accretion).

This is important to determine likely quenching mechanisms. There has previously been a debate about whether major mergers can be involved in quenching massive galaxies (e.g. Ilbert et al. 2013). In this scenario, a major merger would funnel gas into the centre of the galaxy, triggering a starburst, which would rapidly use up the gas and simultaneously cause strong AGN feedback, thereby quenching the galaxy. However, this scenario has been disfavoured (e.g. Grogin et al. 2005; McAlpine et al. 2017; Rodríguez Montero et al. 2019). Overall, while major mergers are likely to cause significant starbursts, they are no longer thought to play a significant role in actually quenching massive galaxies. This is supported by the findings of D’Eugenio et al. (2023) for an early quenched galaxy that could be studied in detail: the quenching process left the dynamics and morphology of the galaxy unaffected (disc-like), implying that it could not have happened via a catastrophic merging event.

According to the simulation predictions, most massive quiescent galaxies at |$z>2\!-\!5$| did not undergo a major merger event recently. This is in qualitative agreement with the compact sizes of observed galaxies, which suggest no major contribution of mergers. Therefore, the quenching of these massive galaxies is most likely related to an internal process. For galaxies of this mass, that is almost definitely going to be related to some form of AGN feedback (e.g. Man & Belli 2018).

5.7 AGN incidence

The duration of active phases in galaxies is debated (Alexander & Hickox 2012; Harrison et al. 2018; Harrison & Ramos Almeida 2024), but time-scales are generally thought to be much shorter than the quenching time-scale of galaxies, of order 1–10 Myr. In fact, both theory and observations suggest that AGN-driven quenching is related to the time-integrated accretion history, as traced by the SMBH mass (e.g. Davies et al. 2020; Bluck et al. 2022, 2023; Brownson et al. 2022; Piotrowska et al. 2022; Baker et al. 2024). Nevertheless, by observing quiescent galaxies at earlier and earlier epochs, we necessarily ‘squeeze’ the timescales of quenching and AGN activity. Specifically, there is less time available over which to spread the time-integrated accretion history, meaning we have a high chance to detect AGN activity. To ensure the most complete census possible, we use a combination of four diagnostics to probe AGN. The first is X-ray detections based on Chandra Legacy data of GOODS-S and GOODS-N fields. For GOODS-N, we use the catalogue and classifications from Xue et al. (2016), whilst for GOODS-S we use the catalogue and classifications from Luo et al. (2017). We cross-match the RA and Dec. to within 1 arcsec. This gives us three matches, all of which are classified as AGN. These are galaxies 199 773, 197 911 (as also identified in Circosta et al. 2019; D’Eugenio et al. 2023; Scholtz et al. 2024), and 83 288. However, non-detection in X-rays does not rule out the presence of an AGN, due to the possibility of substantial obscuration, or even intrinsically weak X-rays emission (e.g. Circosta et al. 2019; Lyu et al. 2024; Juodžbalis et al. 2024a; Mazzolari et al. 2024a, b; Maiolino et al. 2025). This has been ascribed to many different physical processes, but these are beyond the scope of this paper; we refer the reader to Maiolino et al. (2025) for an in-depth summary.

In addition to X-ray selection, we use the classical optical BPT diagram (Baldwin, Phillips & Terlevich 1981) with diagnostic dividing lines from Kauffmann et al. (2003b). We used the flux from the R1000 gratings to separate out H|$\alpha$| and [N ii]. The results are shown in Fig. 12. Note that several spectroscopically confirmed AGN have been reported to occupy the star-forming region of the BPT diagram at redshifts |$z\gtrsim 5$| (e.g. Übler et al. 2023; Maiolino et al. 2024). However, only one of our BPT measurements occupies a region consistent with high-redshift AGN (83 288), and this galaxy is at |$z=2.2$|⁠, below the redshift where the BPT becomes unreliable (Scholtz et al. 2023). We find several AGN candidates which are 011252, 072127, 1171147, 197911, and 199733. Three of these are classified on the basis of the four emission lines, while another two are classified as Low-ionisation Emitting Regions (LIER)/AGN based on their high [N ii]/H|$\alpha$| ratios only.

Two of our sample show broad H|$\alpha$| lines corresponding to a BLR. These are 8777 (previously identified in Carnall et al. 2023b) and 11 252.

We also use the near-infrared (NIR) to help identify AGN. We use the JWST Mid-Infrared Instrument (MIRI) data where available. This identifies 11 252, 197 911, and 199 773 as NIR-selected AGN (all of which are selected as AGN by other diagnostics).

Table 2 provides a summary of our AGN selection criteria. We find we observe AGN activity in 8/18 massive quiescent galaxies in our sample.

5.8 Ly|$\alpha$| emission in a redshift 4.1 quiescent AGN host galaxy

072 127 is a particularly interesting target (Kokorev et al. 2024; Jones et al. 2025). It is extremely compact, with a half-light radius of just 0.3 kpc, comparable to 8777 (Carnall et al. 2023a). It shows clear signatures of being an AGN host (like 8777, as further shown in Ji et al. 2024b), but it displays prominent narrow-line emission, making it a type-2 AGN (unlike 8777, which shows only weak narrow-line emission). Kokorev et al. (2024) recently reported broad-H|$\alpha$| emission, indicating a type-1 AGN, which would further confirm our case for an AGN. However, we note that our analysis does not confirm a strong detection (Juodžbalis, in preparation). In addition to narrow-line optical emission, 72 127 also shows extremely strong Ly|$\alpha$| emission in the prism spectrum (Fig. 1 and Jones et al. 2025), the second such instance in a high-redshift quiescent galaxy after RUBIES-EGS-QG-1 (Urbano Stawinski et al. 2024).

5.9 Morphology and sizes

Another important aspect to explore is the morphology of high-z quiescent galaxies (e.g. see Ferreira et al. 2023; Kartaltepe et al. 2023; Baker et al. 2025a). In the more local to redshift 2 Universe, galaxy sizes and Sérsic indices are important tests of both theory and cosmological models (for a review, see Conselice 2014). This is important to link concepts such as merger rates, growth, and size evolution. As described in Section 3.3, we use PySersic to fit the galaxy images in the F444W filter. We chose this filter because it is the least contaminated by emission lines from the AGN. The downside is that we are most affected by the PSF at this wavelength. To test these effects we also fit the galaxies in the F200W filter and find that the results for the sizes remains consistent.

Fig. 13 shows the effective radius in kpc versus the stellar mass for the 18 massive quiescent galaxies. For comparison, the size–mass relation from der Wel et al. (2014, for early-type galaxies at redshift 2.75) is overplotted. We find that the relation almost completely bisects the galaxies in our sample. Interestingly, we do not appear to see two distinct size–mass relations based on stellar mass (as was seen in a larger photometric sample in Ji et al. 2024a), but this is likely due to the smaller size of our spectroscopically confirmed sample.

We find a large amount of scatter in our size measurements. Despite this, we use orthogonal distance regression to provide a best fit to our sample. We find a shallower best-fitting relation than that of der Wel et al. (2014), but warn that our sample contains higher redshift galaxies and we expect the size–mass relation to have a redshift dependence. Therefore, we caution that the exact shape of our size–mass relation is likely not statistically significant. The relation obtained is |$\rm r_{eff}/kpc= 10^{-0.11\pm 0.08}\times (M_\star /5\times 10^{10}M_\odot)^{0.44\pm 0.21}$| with a scatter of 0.4 dex. We include the effective radii for the quiescent galaxies in Appendix Table B1.

In addition to galaxy sizes, we can also explore the distribution of Sérsic indices for the sample. Quiescent galaxies traditionally show more concentrated light distributions corresponding to Sérsic indices greater than 2. Often these take the form of the traditional de-Vaucouleurs elliptical-style galaxies. We again draw our Sérsic indices from fits to the quiescent galaxies in the F444W filter. We find a range of Sérsic indices from two to eight (eight being the upper bound for the fits). The majority of the galaxies show high Sérsic indices, consistent with them being ellipticals. The Sérsic indices are reported in Appendix Table B1.

6 DISCUSSION

We selected a sample of 18 massive quiescent galaxies from the JADES survey and from PID 2198 (Barrufet et al. 2021) that have both NIRCam photometry and NIRSpec prism spectroscopy, complemented by medium-resolution grating spectroscopy from JADES and PID 1207 for two galaxies. We demonstrated that these galaxies are both quiescent and massive. Our primary goal is to study internal quenching mechanisms; for this reason, we consider only central galaxies, based on the absence of massive neighbours. We reject two lower and one higher mass quiescent satellite galaxies on the basis of their proximity to higher mass galaxies at the same redshift (Section 4.1, and for more details on the satellites, see Appendix  A).

6.1 Number densities

The number density of quiescent central galaxies is one of the most basic tests of galaxy evolution and cosmology. The existence of massive, quiescent galaxies at |$z>3$| is no longer considered a problem for cosmology (e.g. Glazebrook et al. 2017). However, their existence and abundance provide compelling constraints on our models of star formation efficiency (Carnall et al. 2024; Glazebrook et al. 2024; Graaff et al. 2025) and SMBH-driven quenching (Xie et al. 2024). Although a few objects may still be considered outliers (e.g. Carnall et al. 2023b; Glazebrook et al. 2024; Graaff et al. 2025), samples of quiescent galaxies require a revision of our models (Barrufet et al. 2021; Valentino et al. 2023; Carnall et al. 2023a, 2024; Nanayakkara et al. 2024; Park et al. 2024). In this study, we used the superior precision of spectroscopy to study quiescent galaxies down to a stellar mass limit of |$M_\star \gtrsim 10^{10}~\mathrm{M}_\odot $|⁠. In Section 5.1, we compared the number densities of quiescent galaxies in our spectroscopically corrected photometric sample with that found in FLAMINGO, one of the largest box hydrodynamical cosmological simulations to date (Kugel et al. 2023; Schaye et al. 2023). This is a crucial check – many works have previously explored the number densities of high-z quiescent galaxy candidates (as identified photometrically, Valentino et al. 2023; Carnall et al. 2023a; Alberts et al. 2024; Ji et al. 2024a), but we are the first to combine a large, spectroscopically corrected photometric sample with a large-volume hydrodynamic simulation. Our approach removes both possible interlopers and contaminants and enables us to rule out a major role of cosmic variance (e.g. Carnall et al. 2023a). We are also able to account for quiescent galaxies missed photometrically with UVJ colour selection (see Section 4 and Fig. 2).

In addition, the large simulation box allows us to safely probe cosmic variance up to 3|$\sigma$|⁠, as demonstrated in Fig. 5. Our observational analysis reveals a significant excess of massive quiescent galaxies compared to FLAMINGO, both at |$z=2\!-\!3$| (2|$\sigma$|⁠), and even more so at |$z=3\!-\!4.5$| (3|$\sigma$|⁠). This cements the findings of earlier studies (which used photometry and a range of other numerical simulations; e.g. Valentino et al. 2023; Lagos et al. 2025) and gives us important information on the feedback prescriptions of the simulations. This comparison can be seen in more detail in Section 5.2, where we compare our spectroscopically corrected number densities to other observational studies, cosmological simulations and SAMs. We again find that our observed number densities are generally consistent with other observations at high-z, showing an overabundance compared to the cosmological simulations. In particular, our findings suggest that theoretical models require (1) stronger and/or more efficient quenching of massive galaxies, and (2) acting as early as redshifts of 6 and earlier. Energy arguments suggest that this earlier quenching is most likely due to AGN feedback, given that other mechanisms struggle to counteract cosmic gas accretion (e.g. Silk & Rees 1998). We also note the usual caveat that increased resolution would be important to fully probe the parameter space of the typical quenching criteria.

6.2 AGN feedback

That AGN feedback may be stronger and/or more efficient at earlier epochs is suggested not only by the unexpectedly high number density of quiescent galaxies but also by other lines of evidence. In our sample, we find a large fraction of AGN, suggesting continued fuelling of the central SMBH even after hundreds of Myr of quiescence. This is confirmed by several other findings of nebular emission in quiescent galaxies, both at |$z>4$| (e.g. Carnall et al. 2023a; Graaff et al. 2025), in samples near Cosmic Noon (Park et al. 2024; Bugiani et al. 2025), and even at later epochs (Maseda et al. 2021, though these authors argue that low-level star formation could explain their findings). Our work is the first sample study at |$z=2\!-\!5$|⁠, finding an AGN fraction of 50 per cent. A comparison to the AGN incidence in a control sample of star-forming galaxies of the same mass is urgently needed (initial estimates hover around 20 per cent, Scholtz et al. 2023).

We stress that finding AGN activity is not a smoking gun for AGN-driven quenching; simple logic posits that ongoing AGN cannot have caused quenching hundreds of Myr earlier, especially given the short duration of AGN episodes, which is thought to be a few Myr at most (at least at lower redshifts; e.g. Harrison et al. 2018). Evidence from statistical studies at |$0 < z < 2$| finds that quiescence is most closely related to the central stellar velocity dispersion (Bluck et al. 2022, 2023; Baker et al. 2024), which is a proxy for SMBH mass (e.g. Magorrian & Tremaine 1999; Ferrarese & Merritt 2000; Kormendy & Ho 2013), itself a measure of the time-integrated SMBH activity. These studies suggest that it is not any particular AGN episode that quenches a galaxy, but the cumulative effect of all episodes, which heat the galaxy halo and prevent further gas accretion (preventive feedback) (Baker et al. 2024). This finding and its interpretation is supported both by studies of small samples of local galaxies with directly measured SMBH masses, as well as by theoretical predictions (Davies et al. 2020; Bluck et al. 2022; Piotrowska et al. 2022). Recently, evidence in this direction has been reported even by observations of molecular gas (or lack thereof) in a quiescent galaxy at |$z=3$| (Scholtz et al. 2024). Nevertheless, even in this scenario preventive-feedback scenario, we would still expect a higher incidence of AGN at high redshift, simply because at earlier cosmic times, there is less time available to distribute the time-integrated accretion. In addition, maintenance-mode feedback may be required to maintain galaxies quiescent over long periods of time – particularly at early cosmic epochs, when gas density and mass accretion rates were higher. However, we cannot exclude the possibility that AGN-driven quenching was different in earlier epochs (e.g. Xie et al. 2024; Maiolino et al. 2024). Although existing studies based on JWST imaging confirm the trends found at lower redshifts (Bluck et al. 2023), these early studies are still dominated by quiescent galaxies at Cosmic Noon, while quiescent galaxies at earlier epochs may still be too few in number to detect signatures of different trends, as supported e.g. by numerical simulations (Valenzuela & Remus 2024; Remus & Kimmig 2025).

Another independent line of evidence supporting stronger and/or more efficient AGN feedback at high redshift is the discovery of widespread neutral-gas outflows in quiescent galaxies at Cosmic Noon (Davies et al. 2024) and even earlier epochs (Belli et al. 2024; D’Eugenio et al. 2024). These studies confirm the continued fuelling of SMBHs in quiescent galaxies, and underscore the role of ejective feedback – at least in maintaining galaxies as quiescent.

A final possibility is that SMBH growth is more rapid than expected at high redshift. This is supported by the discovery of SMBHs that are ‘overmassive’ relative to the stellar mass of their host galaxy, when compared to local scaling relations (Kocevski et al. 2023; Kokorev et al. 2023; Übler et al. 2023; Mezcua et al. 2024; Maiolino et al. 2024; Juodžbalis et al. 2024b; but see e.g. Sun et al. 2025 or Stone et al. 2024 for a different view). These overmassive black holes are a strong sign that AGN feedback in the early Universe could have been more rapid than what we initially expected, as it outpaces the growth rate of the host galaxy. Alternatively, we could interpret this finding as further evidence supporting stronger/more efficient SMBHs feedback, which may have inhibited star formation in the host galaxy.

Studies of the cold gas reservoir of quiescent galaxies will be crucial to pin down how exactly quenching works (Whitaker et al. 2021; Williams et al. 2021) – particularly at |$z>2\!-\!3$| (Scholtz et al. 2024).

6.3 Star formation histories

Exploring the SFHs of our sample we see a clear divide between galaxies that appear to have formed all their mass within a single short burst and those that appear to have a much more extended star formation period. The ‘burst’ ones, similar to those observed by Carnall et al. (2024), appear to have managed to form hundreds of solar masses worth of stars per year for a short period of time. Something then appears to have rapidly shut off this star formation, and these galaxies have remained dormant ever since. Meanwhile, the more extended SFH galaxies appear to have had a much more gradual decline in star formation, corresponding to a likely much slower acting feedback process.

This kind of divide between fast- and slower-quenching mechanisms has previously been seen for populations of lower redshift quiescent galaxies (e.g. Wu et al. 2018; Belli et al. 2019).

To some extent this can help inform us as to possible quenching mechanisms. In the case of a single AGN outflow being the dominant driver we would expect a very short quenching time-scale and much more of a ‘single-burst’ SFH. Whereas it is possible to envisage the cumulative effects of AGN feedback (e.g. Bluck et al. 2022; Piotrowska et al. 2022; Baker et al. 2024) working to quench a galaxy over an extended period of time. In addition, according to their SFHs, many of our galaxies have remained quenched for extended periods of time, in some cases since z |$\sim$| 6.

This requires a mechanism to prevent further star formation and continue the galaxies’ quiescent phase. Such mechanism cannot be star formation feedback, because even relatively small star-formation episodes would be easily visible due to their low M/L ratios. In contrast, preventative AGN feedback may be able to stop gas accretion and thus explain the protracted quiescent phase. Other quenching mechanisms are possible, e.g. stellar feedback, cosmic rays, environment and more, but these are generally considered less likely quenching causes in these most massive central high-z quiescent galaxies, but have been shown to possibly be important in lower mass high-z galaxies (Gelli et al. 2025).

Our results match nicely with the results of Park et al. (2024) who explored the quenching of massive galaxies at Cosmic Noon with the Bluejay survey. They found evidence for three different types of SFHs: old galaxies quenched at early epochs; galaxies rapidly and recently quenched after a flat/bursty formation history; and galaxies that were quickly and recently quenched after a short and intense burst of star formation (starburst). We appear to find all three scenarios within our sample. We find the early-formed, early-quenched galaxies, which in our case are dominated by the 4000-Å break galaxies (e.g. 11 252 and 72 396), plus galaxies quickly quenched after a longer SFH (corresponding to our extended SFH galaxies, e.g. 199 773 and 27 461) and also the short starbursts (e.g. 80 660, 25 511, or 8777 as found also by Carnall et al. 2023b). However, in addition to the three types found by Park et al. (2024), we also find evidence for galaxies that appear to have quenched more slowly, such as 83 288 and 200 733.

The formation and quenching epochs within our high-z galaxy sample are generally at much higher redshifts consistent with results found for individual galaxies and much smaller samples (e.g. Carnall et al. 2023b; Glazebrook et al. 2024; Weibel et al. 2024a; Graaff et al. 2025). Part of the reason is due to selection; our sSFR cut removes more recently quenched galaxies (i.e. those having quenched within the 100 Myr prior to observation). However, despite this ‘adverse’ selection, we still find many more quiescent galaxies at these redshifts than would have been expected in the pre-JWST era (Section 6.1), in agreement with earlier studies (Valentino et al. 2023; Carnall et al. 2023a).

As part of our analysis we did explore whether we saw any correlation between formation and quenching time and other galaxy properties such as stellar mass for our 18 quiescent galaxies. This can probe important galaxy effects such as’downsizing’ (Cowie et al. 1996; Thomas et al. 2005; Ilbert et al. 2013). We see a tentative trend towards more massive galaxies having an older formation and quenching time, but due to the small nature of our sample, we find that it is not statistically significant. Larger samples of spectroscopically confirmed high-z massive quiescent galaxies will be required to properly test these trends.

We find no evidence for rejuvenating galaxies (Witten et al. 2025), or galaxies that appear ‘mini’ or ‘rapidly’ quenched (Looser et al. 2023; Dome et al. 2024; Looser et al. 2024; Baker et al. 2025b), although again we caution that this is due to our selection criteria. In particular, while some models predict that quiescent galaxies at high redshift are likely to rejuvenate (Remus & Kimmig 2025), our data may not be sensitive to bimodal SFHs if both peaks are older than 100 Myr (as required by our selection criteria). Investigating recently rejuvenated galaxies is easier when including currently star-forming galaxies, which increases the baseline in logarithmic lookback time over which the two peaks of SFR can be separated. As for mini or rapidly quenched galaxies, medium-band photometric studies find that these systems may be abundant in GOODS-S (Trussler et al. 2025). However, they are by definition rapidly and recently quenched galaxies, with quenching times of a few tens of Myr (e.g. Looser et al. 2024; Baker et al. 2025b) – which we specifically select against. In addition, mini-quenched galaxies are expected to have a much lower mass than those we are probing and to be more common at higher redshifts than what we consider here (Ceverino, Klessen & Glover 2018; Ceverino et al. 2021; Lovell et al. 2023; Dome et al. 2024). For these reasons, we find no tension between our lack of mini-quenched systems and previous studies (e.g. Trussler et al. 2025).

Interestingly, four galaxies within our sample exhibit 4000-Å breaks, an unequivocal indicator of an aged stellar population. These breaks are distinct from the Balmer break, which typically describes stellar populations that are a few hundreds of Myr old. In contrast, the 4000-Å break typically appears in populations around 1 Gyr old and older, depending on metallicity.

These 4000-Å break galaxies have been seen before in high-z studies, most notably in the case of ZF-UDS-7329 (Glazebrook et al. 2024). This galaxy was observed at a redshift of 3.205, when the Universe was around 2 Gyr old. The key conundrum with these breaks is that they imply a much earlier formation time, as they are only seen in old stellar populations, which provides strong constraints to galaxy formation. Indeed, Glazebrook et al. (2024) show that there is a lack of massive enough dark-matter haloes to host their galaxy at the redshift at which it should have built up its measured stellar mass.

The galaxies in our sample that show 4000-Å breaks appear to be less extreme than the Glazebrook et al. (2024) example. They are all at lower redshifts and less massive. However, the fact that we see four of them suggests that they may be much more prevalent than previously thought. Their formation epochs of |$z>8$|⁠, coupled with early quenching around redshift 5–8, require rapid quenching. They also require no recent star formation, which would otherwise erase the 4000-Å break (similar to what is described for Balmer breaks; Witten et al. 2025). Fascinatingly, even in the 4000-Å break subset, two of them turn out to be AGN hosts (in agreement with the overall sample fraction of 50 per cent; Section 6.2). This suggests that there is still some form of gas reservoir/fuel remaining, not only in recently quenched galaxies, but even in these maximally old galaxies.

6.4 In-situ versus ex-situ star formation

Using the FLAMINGO simulation, we were able to explore the mass build-up of massive quiescent galaxies, finding that primarily the mass is formed in situ (Section 5.6), with only a small fraction of ex-situ stars. The simulations predict that these galaxies typically experience one or two major dry mergers, or only 0.6 Gyr⁻¹. They can still undergo minor mergers that may contribute more than major mergers to the mass budget. However, the numerical resolution of FLAMINGO does not allow us to probe down to minor mergers. This result, in combination with the high AGN incidence (Section 6.2), presents an interesting scenario. 50 per cent of our massive quiescent galaxies host an AGN (and the true number is likely to be higher). This requires the presence of gas to fuel the AGN, but this gas cannot have come from the classic gas-rich, major-merger scenario, which would efficiently channel gas into the galaxy’s centre. This leaves minor mergers or halo cooling as sources of gas.

However, as previously mentioned, there remains the possibility that we cannot accurately pick up wet mergers due to star formation resulting from the gas of the wet merger. This would then be counted as in-situ star formation, possibly washing out any sudden increase in the stellar mass budget. Alternatively, AGN feedback is ineffective in expelling all of the gas from the galaxies. This would leave a significant gas reservoir that, after spending some time at large galactocentric distances, could return to continue to fuel the AGN.

6.5 Limitations of the simulations

Although we exploited the large box size of FLAMINGO simulations as a crucial advantage to explore the cosmic variance, it comes at several expenses as well, including the resolution. Increasing a numerical resolution normally leads to a lower quenched fraction, because of more star forming gas being resolved, and assignment of non-zero SFR to galaxies whose SFR fall below a poorer resolution. Therefore, given the direction of the current discrepancy, namely fewer quenched objects predicted by the simulations than observations, our simulation results for the number densities can be considered rather as a conservative limit to the tension in principle. For example, the formation of a single star particle in 100 Myr corresponds to a sSFR = 1.3|$\times 10^8$|  |$\rm M_\odot$|/|$100\times 10^8 \rm yr$|/|$\rm 2.5\times 10^{10}\,\,M_\odot$| = |$5.2\times 10^{-11}\,\,\rm yr^{-1}$| at the observed lower mass limit. Our sSFR quenching criterion corresponds to a sSFR |$\rm \le 0.2/3\times 10^9\,\,yr^{-1}$| = |$\rm 6.7\times 10^{-11}\,\,yr^{-1}$| at |$z=2$|⁠, which is clearly not well resolved at the lower mass limit.

However, there are a few specifics of FLAMINGO that complicate such interpretation and the impact of resolution. First, the simulations take so-called weak convergence scheme where the values of the model parameters are adjusted when changing the numerical resolution. This means that the simulations were tuned to best match the observational constraints at each resolution independently. Although this is generally good, because one does not need to care about an absolute numerical convergence, it complicates assessing the significance and even the direction of resolution effect. Namely, it is practically impossible to know whether more or less quenched galaxies will be found when the resolution increases. Second, the simulations fit to the SMF at |$z\, {\simeq }\, 0$| only up to |$\rm M_\ast \, {=}\, 10^{11.5}\, {\rm M}_\odot$|⁠, to calibrate the model parameters. While this is due to the great uncertainties of the SMF at the massive end, this may impact the reliability of the predictions assuming a large fraction of our samples are possibly progenitors of those massive galaxies not used by the simulations as constraints.

There is another limitation of using simulations, in general, for studying rare, high-z objects. Many of the state-of-the-art hydrodynamical cosmological simulations use |$z\, {\simeq }\, 0$| observations of normal populations of galaxies as primary constraints for the calibration of their models, as is also the case for FLAMINGO. This complicates assessing the accuracy and reliability of their predictions on rare, extreme galaxies at high redshifts, since there could be a wide range of evolutionary paths that can lead to very similar low-redshift scaling relations, hence degeneracy. We refer the reader to Lim et al. (2021) for a more extensive discussion of the limitations involved in comparisons of high-z extreme objects between observations and models in general.

7 CONCLUSIONS

We explored spectroscopically confirmed massive quiescent galaxies from the JADES survey to investigate their SFHs, quenching time-scales, and AGN incidence. We also explored a representative sample from the FLAMINGO simulations to compare number densities and mass build-up. Our key results are the following.

• We find 18 massive spectroscopically confirmed quiescent galaxies, five of which are above redshift 4, some of which are clearly quiescent but do not fulfill the traditional UVJ cut (but are still selected via a sSFR criterion).

• We observe more spectroscopically confirmed massive quiescent galaxies than predicted from the FLAMINGO simulations, even more than allowed by the 3|$\sigma$| cosmic variance (especially at |$z=3\!-\!4.5$| and |$4\!-\!5$|⁠). This has important implications for theoretical modelling, in particular for the feedback prescriptions used in the simulations.

• We find both short and long SFHs suggesting the prevalence of both short and long time-scale quenching mechanisms.

• We find four galaxies with clear evidence of 4000-Å breaks, a signifier of an old evolved stellar population. These galaxies have high formation and quenching redshifts and are particularly interesting objects.

• We find that most of our massive quiescent galaxies formed above redshift 5, with some of these likely to have formed above redshift 6. The majority of our sample quenched at redshifts 4–6 with a few galaxies having quenched even earlier.

• Using the FLAMINGO simulation, we explore the stellar mass build-up of massive quiescent galaxies similar to our samples. We find that the mass of these galaxies is almost entirely built up in-situ with ex-situ mass (via mergers) contributing little to the total mass budget. This is explained by these galaxies on average appearing to have not undergone a major merger event recently.

• Using a combination of diagnostics, AGN activity is found within 8/18 galaxies in our observed sample. This means that these galaxies have had a source of fuel for the AGN so are likely to contain some non-negligible gas reservoirs.

• The sizes and morphologies of the observed quiescent galaxies are explored. Every quiescent galaxy in our sample appears to be some form of spheroid with a Sérsic index above 2. Many appear to have high-Sérsic indices above 5, showing a compact central light distribution. In many cases, it is likely to trace the AGN component.

ACKNOWLEDGEMENTS

WMB, FDE, RM, TJL, and JS acknowledge support by the Science and Technology Facilities Council (STFC), by the ERC through advanced grant 695671 ‘QUENCH’, and by the UKRI Frontier Research grant RISEandFALL. This work was supported by a research grant (VIL54489) from Villum Fonden. RM also acknowledges funding from a research professorship from the Royal Society. SA acknowledges grant PID2021-127718NB-I00 funded by the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI/ 10.13039/501100011033) AJB acknowledges funding from the ‘FirstGalaxies’ advanced grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 789056). SC acknowledges support by European Union’s HE ERC Starting grant no. 101040227 – WINGS. XJ acknowledges support by JWST/NIRCam contract to the University of Arizona, NAS5-02015, BER acknowledges support from the NIRCam Science Team contract to the University of Arizona, NAS5-02015, and JWST Program 3215. The authors acknowledge use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315. HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship. The research of CCW is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. CNAW acknowledges support from JWST/NIRCam contract to the University of Arizona, NAS5-02015. JZ acknowledges support from JWST/NIRCam contract to the University of Arizona, NAS5-02015. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1, and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. This work is based [in part] on observations made with the NASA/ESA/CSA JWST. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–03127 for JWST. These observations are associated with programs nos 1180, 1181, 1210, 1286, 1895, 1963, and 3215. The data are from the JADES survey with the associated MAST, doi:10.17909/8tdj-8n28.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author. The JADES data are publicly available at https://archive.stsci.edu/hlsp/jades/.

Footnotes

REFERENCES

APPENDIX A: HOW DEEP CAN WE GO? BRIEF EXPLORATION OF QUIESCENT SATELLITE GALAXIES

Although the goal of this paper is to explore the massive quiescent galaxy population within GOODS-S and GOODS-N, in the process of these explorations, we also uncover two lower mass (likely satellite) quiescent galaxies, 10010639 and 7300 in addition to a higher mass satellite galaxy 10026167. These satellite galaxies are shown as RGB images in Fig. A1. Galaxy 7300 has previously been identified and studied in Sandles et al. (2023), so it is reassuring to rediscover it. We can see it is compact and faint in the imaging. Galaxy 10010639 is as of yet unpublished and appears to also be low-mass. We can clearly see a large central galaxy to the south-east in the imaging. However, 10026167 appears to be a higher mass quiescent satellite galaxy with a stellar mass of |$M_\star =10^{10.4}\,\mathrm{M}_\odot $|⁠. In order to check that this is not a central, we fit the photometry of the massive dusty star-forming neighbour 175041 (upper left of the RGB image) in order to obtain a stellar mass estimate. We find 175041 has a stellar mass of |$M_\star =10^{11.3}\,\mathrm{M}_\odot $| so is clearly a massive central. It also shows signs of considerable star-formation activity with an |$\rm SFR_{10}=22^{+8}_{-6}\,M_\odot\,yr^{-1}$|⁠. This means that we can be confident that 10026167 is in fact a quiescent satellite of this central galaxy (that is itself not quiescent). Exact explorations of these three satellite galaxies’ environment is beyond the scope of this project but will be followed up in the future.

Fig. A2 shows the spectra + photometry and SFHs for the two low-mass and one higher mass quiescent satellite galaxies. We can clearly see from the SFHs that all three galaxies are quiescent with no signs of recent star formation. The three appear to have gone through short bursts of star formation lasting roughly 300 Myr before quenching. We find a stellar mass of |$10^{9.4}\,\mathrm{M}_\odot $| for 10010639 and |$10^{8.6}\,\mathrm{M}_\odot $| for 7300, whilst as previously mentioned 10026167 has a stellar mass of |$10^{4}\,\mathrm{M}_\odot $|⁠. These galaxies provide a fascinating avenue for further follow up research to better understand possible environmental quenching.

APPENDIX B: TABLES

Table B1 is a table showing the NIRSpec ID (JADES in the case of the JADES spectra and DJA in the case of the four Barrufet et al. 2021 spectra), right ascension and declination, and the best-fitting properties including redshift, stellar mass, formation redshift (⁠|$z_{50}$|⁠), half-light radius, Sérsic index, Balmer break strength, and D|$_n$|4000 index.

We measure the Balmer break strength in |$F_\nu$| following the procedure in Roberts-Borsani et al. (2024) and Witten et al. (2025). Briefly, this is as follows: we use a window between 3500–3630 Å (⁠|$F_\nu (3565\,\rm{\mathring{\rm A}})$|⁠) and a window between 4160–4290 Å (⁠|$F_\nu (4225\rm{\mathring{\rm A}})$|⁠) in order to select regions that do not overlap with strong emission lines. We then take the ratio between the two to give the Balmer break strength such that

The D_n4000 index is measured as defined in Balogh et al. (1999).

APPENDIX C: EXTRA SPECTRA

Figs C1 and C2 contain the observed spectra and photometry and best-fitting model spectra and photometry for the remaining nine galaxies in our sample. The other nine galaxies’ spectra and photometry appeared in Figs 1 and 10. Once again, the observed spectra and errors are in grey, and the observed photometry is given by the yellow points. The best-fitting model spectra is in red, and the best-fitting model photometry is given by the black squares. The lower panel shows the |$\chi$| values of the fit to the spectra.

APPENDIX D: EXTRA STAR FORMATION HISTORIES

Figs D1 and D2 show the SFHs for the remaining 12 massive quiescent galaxies. The vertical dashed lines correspond to the formation time |$\rm t_{50}$| and the quenching time |$\rm t_{90}$|⁠.

APPENDIX E: FITTING BOTH PHOTOMETRY AND SPECTRA WHEN OFFSET

Fig. E1 shows the NIRSpec prism spectrum (grey) and NIRCam photometry (yellow points) for the massive quiescent galaxy 171 147. This galaxy is significantly larger than the NIRSpec slit, as can be seen in Fig. 3. Therefore, it provides a helpful example of how we fit both photometric and spectroscopic data at the same time. We want to include information from the spectrum, whilst also being aware that this traces a subregion of the galaxy, whereas the photometry traces the overall galaxy as a whole. The red spectrum in Fig. E1 corresponds to our best-fitting model of the spectra and the photometry. We see that our best-fitting model accurately traces the shape and flux of the photometry (i.e. the galaxy as a whole) but is able to also take into account the superior information from the spectrum on quantities such as redshift, Balmer break shape, emission lines, and more. Clearly, this methodology is not a substitute for missing data; there could be scenarios where the spectrum differs substantially from the photometry for physical reasons (e.g. the spectrum captures only the quiescent bulge of an early spiral, while the star-forming arms are completely outside the NIRSpec shutter). However, the compact nature and smooth morphology of most of our targets suggest that such effects are not a concern in this work.

APPENDIX F: FITTING THE GRATING SPECTRA

An interesting test is to explore what effects fitting the grating spectra in addition to the prism spectrum has on the SFH. This enables us to explore the effects of higher resolution data which should be better able to fit absorption features within the quiescent galaxies spectrum. We test this by using galaxy 11 252 for which we have gratings spectra. Fig. F1 shows the best-fitting spectra to the prism, grating, and photometry and the difference between the SFHs for the grating + prism and prism only fits. We see that both SFHs are consistent to within 1|$\sigma$|⁠. We do obtain differences in the formation and quenching times, with them being roughly 30 per cent older with the inclusion of the grating data. This highlights the effects of additional grating data, particularly in the case of AGN with strong emission lines.

APPENDIX G: PySersic FITS

Fig. G1 shows the data, model, and residual maps for the F200W morphological fits of three quiescent galaxies in our sample. The galaxies (and additional sources) appear to be well modelled by the single Sérsic profile fits without leaving significant residuals.