# A surprising abundance of massive quiescent galaxies at $\mathbf{3<z<5}$ in
the first data from JWST CEERS

A. C. Carnall1, D. J. McLeod1, R. J. McLure1, J. S. Dunlop1, R. Begley1, F.
Cullen1, C. T. Donnan1, M. L. Hamadouche1, S. M. Jewell1, E. W. Jones1, C. L.
Pollock1, V. Wild2
1 SUPA, Institute for Astronomy, University of Edinburgh, Royal Observatory,
Edinburgh EH9 3HJ, UK
2 SUPA††footnotemark: , School of Physics and Astronomy, University of St
Andrews, North Haugh, St Andrews KY16 9SS, UK E-mail<EMAIL_ADDRESS>Universities Physics Alliance

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We report a robust sample of 10 massive quiescent galaxies at redshift, $z>3$,
selected using the first data from the JWST CEERS programme. Three of these
galaxies are at $4<z<5$, constituting the best evidence to date for quiescent
galaxies significantly before $z=4$. These extreme galaxies have stellar
masses in the range log${}_{10}(M_{*}/$M${}_{\odot})=10.1-11.1$, and formed
the bulk of their mass around $z\simeq 10$, with two objects having star-
formation histories that suggest they had already reached
log${}_{10}(M_{*}/$M${}_{\odot})>10$ by $z\gtrsim 8$. We report number
densities for our sample, demonstrating that, based on the small area of JWST
imaging so far available, previous work appears to have underestimated the
number of quiescent galaxies at $3<z<4$ by a factor of $3-5$, due to a lack of
ultra-deep imaging data at $\lambda>2\,\mu$m. This result deepens the existing
tension between observations and theoretical models, which already struggle to
reproduce previous estimates of $z>3$ quiescent galaxy number densities.
Upcoming wider-area JWST imaging surveys will provide larger samples of such
galaxies and more-robust number densities, as well as providing opportunities
to search for quiescent galaxies at $z>5$. The galaxies we report are
excellent potential targets for JWST NIRSpec spectroscopy, which will be
required to understand in detail their physical properties, providing deeper
insights into the processes responsible for forming massive galaxies and
quenching star formation during the first billion years.

###### keywords:

galaxies: evolution – galaxies: star formation – methods: statistical

††pubyear: 2022††pagerange: A surprising abundance of massive quiescent
galaxies at $\mathbf{3<z<5}$ in the first data from JWST CEERS–8

##  1 Introduction

Two of the most important outstanding questions in galaxy evolution are: when
did the first galaxies begin to form stars, and when did the first galaxies
quench their star-formation activity? During the short time since the first
data from the James Webb Space Telescope (JWST) were released, remarkable
progress has been made towards addressing the first of these questions. We now
have good evidence that log${}_{10}(M_{*}/$M${}_{\odot})\simeq 8-9$ galaxies
were already in place by $z\simeq 17$, less than 250 Myr after the Big Bang,
with preliminary evidence mounting that such objects are more numerous than
expected (e.g. Naidu et al. 2022; Castellano et al. 2022; Donnan et al. 2023;
Finkelstein et al. 2022). We have also uncovered the unexpectedly rapid growth
of these early seeds into massive galaxies with
log${}_{10}(M_{*}/$M${}_{\odot})\simeq 10-11$ during the latter half of the
first billion years, from $6<z<10$ (Labbe et al., 2022).

This extremely rapid assembly of the first massive galaxies is critically
important for our understanding of quenching. If
log${}_{10}(M_{*}/$M${}_{\odot})\simeq 11$ galaxies already exist by $6<z<10$,
these must equally rapidly quench, and remain quenched, to avoid becoming too
massive to be accommodated by the lower-redshift galaxy stellar-mass function
(e.g. McLeod et al. 2021). This suggests massive quiescent galaxies at least
as early as $z\simeq 6$.

Currently, the earliest spectroscopically confirmed massive quiescent galaxies
are half a billion years later, when the Universe was $\simeq 1.5$ Gyr old at
$z\simeq 4$ (e.g. Glazebrook et al. 2017; Valentino et al. 2020; Forrest et
al. 2020), and indeed it has proven extremely challenging to identify even
robust photometric candidates at $z>4$ (e.g. Merlin et al. 2018, 2019; Carnall
et al. 2020; Esdaile et al. 2021; Santini et al. 2021; Stevans et al. 2021;
Marsan et al. 2022). It is currently unclear whether this is due to an almost
total lack of quiescent galaxies at earlier times, or due to a lack of ultra-
deep, high-resolution imaging at $\lambda>2\,\mu$m with which to constrain the
Balmer break at these redshifts.

The number density and passive fraction of high-redshift massive galaxies are,
however, key constraints on galaxy formation models, with current simulations
unable to reproduce the observed number density of quiescent galaxies at
$3<z<4$ (e.g. Schreiber et al. 2018; Cecchi et al. 2019; Girelli et al. 2019).
This implies that key physics, capable of giving rise to extremely rapid
quenching events, is still missing from these simulations. In this context,
more-robust constraints on the number density of massive quiescent galaxies at
$3<z<4$, and confirmation of whether any such objects exist at $z>4$, are key
to our understanding of galaxy formation.

Figure 1: Spectral energy distributions and cutout images for our three robust $z>4$ quiescent galaxies. Our 10-band photometric data from HST ACS and JWST NIRCam are shown in blue and gold respectively. The posterior median Bagpipes models are overlaid in red. Posterior distributions for the redshifts and sSFRs of these galaxies are shown to the right of the main panels. The dashed vertical lines in the sSFR panels show the sSFR threshold for inclusion in our quiescent sample at the redshift of each object (see Section 3.2). The inset RGB cutouts are composed of the F444W, F200W and F150W images respectively. Table 1: The 9 free parameters of the Bagpipes model we fit to our photometric data (see Section 3), along with their associated prior distributions. The upper limit on $\tau$, $t_{\mathrm{obs}}$, is the age of the Universe as a function of redshift. Logarithmic priors are all applied in base ten. For the Gaussian prior on $\delta$, the mean is $\mu$ and the standard deviation is $\sigma$. Component | Parameter | Symbol / Unit | Range | Prior | Hyper-parameters
---|---|---|---|---|---
General | Redshift | $z$ | (0, 20) | Uniform |  |
| Total stellar mass formed | $M_{*}\ /\ \mathrm{M_{\odot}}$ | (1, $10^{13}$) | Logarithmic |  |
| Stellar and gas-phase metallicities | $Z\ /\ \mathrm{Z_{\odot}}$ | (0.2, 2.5) | Logarithmic |  |
Star-formation history | Double-power-law falling slope | $\alpha$ | (0.01, 1000) | Logarithmic |  |
| Double-power-law rising slope | $\beta$ | (0.01, 1000) | Logarithmic |  |
| Double power law turnover time |  $\tau$ / Gyr | (0.1, $t_{\mathrm{obs}}$) | Uniform |  |
Dust attenuation |  $V-$band attenuation |  $A_{V}$ / mag | (0, 8) | Uniform |  |
| Deviation from Calzetti et al. (2000) slope | $\delta$ | ($-0.3$, 0.3) | Gaussian | $\mu=0$ |  $\sigma$ = 0.1
| Strength of 2175Å bump | $B$ | (0, 5) | Uniform |  |

By providing unprecedentedly deep infrared imaging at $\lambda>2\,\mu$m, JWST
opens up a unique opportunity to address this issue, promising the ability to
robustly select representative samples of massive galaxies as far back as just
a few hundred Myr after the Big Bang. In addition, its extremely high angular
resolution (e.g. Suess et al. 2022) and wide-ranging spectroscopic
capabilities (e.g. Carnall et al. 2023b) hold much promise for extending
detailed studies of quiescent galaxy physical properties, such as star-
formation histories (SFHs), stellar metallicities and sizes, back to the first
billion years. This endeavour has previously proven extremely challenging,
even at cosmic noon (e.g. Wu et al. 2018; Belli et al. 2019; Carnall et al.
2019a, 2022; Beverage et al. 2021; Hamadouche et al. 2022).

In this paper, we use extremely deep $\lambda=1-5\,\mu$m NIRCam imaging from
the JWST Cosmic Evolution Early Release Science
(CEERS111https://ceers.github.io) programme (Finkelstein et al. in prep.) to
search for massive galaxies at $z>3$. In particular, we focus on constraining
the number density of galaxies that have already quenched their star-formation
activity at this early time. We also make a first attempt at measuring the
SFHs of these galaxies, despite the significant challenge of measuring these
from photometric data alone, in order to link them with the extreme population
of star-forming galaxies currently being uncovered during the first billion
years. The spectroscopic capabilities of JWST mean that spectroscopic
redshifts, SFR measurements, and even spectroscopic SFH and stellar
metallicity determinations are a realistic prospect for these galaxies on a
short timescale (e.g., Carnall et al. 2023a).

The structure of this paper is as follows. In Section 2 we introduce the CEERS
and ancillary datasets used in this work. In Section 3 we discuss our spectral
energy distribution (SED) fitting methodology and sample selection. In Section
4 we present our results, including the discovery of three robustly identified
massive quiescent galaxies at $4<z<5$. We discuss our results in Section 5 and
present our conclusions in Section 6. All magnitudes are quoted in the AB
system. For cosmological calculations, we adopt $\Omega_{M}=0.3$,
$\Omega_{\Lambda}=0.7$ and $H_{0}$ = 70 $\mathrm{km\ s^{-1}\ Mpc^{-1}}$. We
assume a Kroupa (2001) initial mass function, and assume the Solar abundances
of Asplund et al. (2009), such that $\mathrm{Z_{\odot}}=0.0142$.

Figure 2: Comparison of pre-JWST and new JWST CEERS fit results for object
101962, the only robust $z>4$ quiescent galaxy in our sample to be included in
the CANDELS photometric catalogue of Stefanon et al. (2017). The best fit to
the new JWST NIRCam data (gold) is shown in red, as in the middle panel of
Fig. 1. The best fit to the CANDELS data (blue) is shown in green. The corner
plot to the left shows constraints on the stellar mass, redshift and sSFR of
this galaxy from both datasets. It can be seen that previous data were unable
to constrain these parameters, with an extremely large redshift uncertainty,
$z=1.34^{+1.94}_{-0.39}$. This is largely due to the low SNR of previous data,
in particular around the Balmer break, and a lack of data at $\lambda\simeq
2.3-3.0\mu$m between $K_{s}$ and IRAC Channel 1.

##  2 Data

The primary dataset for this work is comprised of the first observations made
as part of the CEERS survey (Finkelstein et al. in prep) in the CANDELS
Extended Groth Strip (EGS) field. We use data from the first 4 pointings
observed in late June 2022. Imaging is available in 7 NIRCam filters: F115W,
F150W, F200W, F277W, F356W, F410M and F444W, with integration times of 2635
seconds per filter, except F115W where the exposure time is doubled. The
currently available CEERS data amounts to a total effective area of $\simeq
30$ square arcmin (e.g. Donnan et al. 2023).

In addition to the NIRCam imaging, we also make use of Hubble Space Telescope
(HST) ACS data in the F435W, F606W and F814W bands. We use the v1.9 EGS
mosaics produced by the CEERS team (Koekemoer et al., 2011), which have a
pixel scale of $0.03^{\prime\prime}$.

We perform our own custom reduction of the NIRCam data, beginning with the
Level 1 data products, using the PRIMER Enhanced NIRCam Image Processing
Library (PENCIL), a custom version of the JWST pipeline (v1.6.2). We align and
stack the individual reduced images using the Scamp and Swarp codes (Bertin,
2006, 2010). Our final mosaic images have a pixel scale of
$0.03^{\prime\prime}$ in all bands. We PSF-homogenise all data to the F444W
band using empirical PSFs derived from stacks of bright stars in our mosaic
images. We make use of the CRDS_CTX = jwst_0942.pmap version of the JWST
calibration files, released on 28th July 2022. We then apply the empirically
derived, module-dependent calibration correction factors described in Appendix
C of Donnan et al. (2023). These calibrations are in good agreement with
empirical calibrations derived by other teams222e.g.
https://github.com/gbrammer/grizli/pull/107.

We employ the SExtractor code (Bertin & Arnouts, 1996) to measure object
photometry. We run SExtractor in dual image mode, with the F200W mosaic used
as the detection image in all cases, as this is the wavelength range in which
the SEDs of $z\simeq 3-5$ quiescent galaxies peak. Fluxes are measured within
$0.5^{\prime\prime}$ (16 pixel) diameter circular apertures. We then aperture-
correct by scaling to measurements of FLUX_AUTO in the F200W band. This size
of aperture is more than sufficient for our purposes, given that massive
quiescent galaxies at high redshift are known to be both extremely compact and
centrally concentrated (e.g. van der Wel et al. 2014).

We cut our catalogue at F200W = 26.5, at which magnitude objects have a
typical F200W SNR $\simeq 20$. This is necessary to ensure objects are
reliably detected with high SNR in the other relevant bands, in particular to
obtain strong constraints on the Balmer break strength. This results in a
parent sample of 10542 objects.

We measure uncertainties for each object as the standard deviation of fluxes
measured in the closest $\simeq 100$ blank sky apertures (McLeod et al.,
2016), whilst masking out nearby objects, using the robust median absolute
deviation (MAD) estimator. We check the resulting photometry in the F606W and
F814W bands by cross-matching with the CANDELS catalogue produced by Stefanon
et al. (2017), finding good agreement (e.g., median offset of 0.07 magnitudes
in F814W, with a scatter of 0.14 magnitudes).

Figure 3: The rest-frame UVJ colour diagram, showing our magnitude-selected
sample of massive galaxies (F200W < 26.5) in integer redshift bins spanning
$3<z<6$. Points are coloured by sSFR, with star-forming galaxies denoted by
circles without a border. Objects in our robust quiescent sub-sample are shown
with black-bordered squares, whereas objects in our quiescent sample that do
not meet our robust criteria (see Section 3.2) are shown as black-bordered
circles. SEDs for the three robust quiescent galaxies we identify at $z>4$ are
shown in Fig. 1, whereas the 7 robust quiescent galaxies at $3<z<4$ are shown
in Fig. 7.

##  3 Method

###  3.1 Spectral energy distribution fitting

The SED-fitting analysis in this work makes use of the Bagpipes spectral
fitting code (Carnall et al., 2018), and is based on the method used to search
for $2<z<5$ massive quiescent galaxies in CANDELS UDS and GOODS South in
Carnall et al. (2020).

The model we fit to our photometric data makes use of the Bruzual & Charlot
(2003) stellar population models, in particular the 2016 updated version
(Chevallard & Charlot, 2016) using the MILES stellar spectral library
(Sánchez-Blázquez et al., 2006; Falcón-Barroso et al., 2011) and the updated
stellar evolutionary tracks of Bressan et al. (2012) and Marigo et al. (2013).

Nebular line and continuum emission are included in our model using an
approach based on the Cloudy photoionization code, outlined in section 3 of
Carnall et al. (2018), following Byler et al. (2017). We assume an ionization
parameter, $U=10^{-3}$, and a lifetime for stellar birth clouds of 10 Myr.

Dust attenuation is included using the model of Salim et al. (2018), which has
a variable slope, parameterised with a power-law deviation, $\delta$, from the
Calzetti et al. (2000) model. We allow the $V-$band attenuation, $A_{V}$, to
vary from $0-8$ magnitudes. We further assume that light from stars still
enclosed in stellar birth clouds and resulting nebular emission is attenuated
by twice the $A_{V}$ experienced by older stars within the wider interstellar
medium (ISM) of the galaxy (e.g. Charlot & Fall 2000).

We assume a double-power-law SFH model, as introduced in Carnall et al. (2018,
2019b), which has been shown to reproduce well the SFHs of massive quiescent
galaxies in the Mufasa simulation (Davé et al., 2016). The stellar and nebular
metallicities of galaxies are assumed to be identical, and are varied with a
uniform prior in logarithmic space from $-0.7<$
log${}_{10}(Z/$Z${}_{\odot})<0.4$. Intergalactic medium absorption is included
using the model of Inoue et al. (2014). We vary redshift in our model with a
uniform prior over the redshift range $z=0-20$. A full list of the 9 free
parameters of our model and their associated prior distributions is given in
Table 1. We fit our Bagpipes model to the data using the MultiNest nested
sampling algorithm (Skilling, 2006; Feroz et al., 2019), accessed via the
PyMultiNest interface (Buchner et al., 2014).

###  3.2 Selection of massive quiescent galaxies

We begin our selection process by requiring that objects have a posterior
median redshift greater than $z=3$. We also require that 97.5 per cent of the
redshift posterior for each object lies above $z=2.75$. This is in order to
exclude objects with significant secondary low-redshift solutions, whilst
retaining objects with narrow redshift posteriors that extend marginally below
$z=3$.

The sample is then cleaned by visual inspection of all 10 photometric bands,
as well as the fitted Bagpipes SEDs. We exclude objects that fall close to the
edges of the NIRCam detector, objects for which coverage is only available in
some bands, and various kinds of NIRCam detector artefacts (Rigby et al.,
2022). We further exclude objects that are visible in short-wavelength
imaging, below the position of the Lyman break at the Bagpipes fitted redshift
(which strongly implies the fitted redshift is incorrect). At the end of this
process, we have a total of 421 galaxies at $z>3$.

To separate star-forming and quiescent galaxies, we use a time-dependent cut
in specific star-formation rate (sSFR), as has been widely applied in the
literature (e.g. Gallazzi et al. 2014; Pacifici et al. 2016). We define
quiescent galaxies as those that have

$\mathrm{sSFR}<\dfrac{0.2}{t_{\mathrm{obs}}}$ (1)

where $t_{\rm obs}$ is the age of the Universe at the redshift of the galaxy.
This threshold is broadly equivalent to a selection in rest-frame UVJ colour
space of $U\ -\ V\ >\ 0.88\ \times\ (V\ -\ J)\ +\ 0.69$ (Carnall et al., 2018,
2019b) at all redshifts, which is the $z<0.5$ quiescent galaxy selection
criterion introduced by Williams et al. (2009).

Our full passive sample is defined as those for which the 50th percentile of
the fitted Bagpipes sSFR posterior distribution falls below this threshold.
These objects are robustly placed at $z>3$ by our fitting, and are more likely
to be quiescent than star forming. However, we cannot confidently exclude
star-forming solutions in all cases. Following Carnall et al. (2020), we then
further define a “robust” quiescent sub-sample, for which 97.5 per cent of the
sSFR posterior is required to fall below the threshold in Equation 1. For
these robust objects, we exclude both low-redshift and star-forming solutions
with high confidence.

By this process we identify a total of 15 objects, 10 of which satisfy our
robust selection criteria. We visually inspect the available $Spitzer$-MIPS
24$\mu$m data (Dickinson & FIDEL Team, 2007) for these objects to check for
anomalously strong detections that would clearly indicate these are lower-
redshift dusty star-forming galaxies, however we do not find any obvious
bright detections.

Table 2: Properties of the 15 massive quiescent galaxies at $z>3$ we identify in this work. We define $t_{\mathrm{form}}$ as the age of the Universe corresponding to the mass-weighted age of each galaxy, and $z_{\mathrm{form}}$ as the corresponding redshift. For objects labelled as “robust”, we can confidently exclude star-forming solutions, whereas for non-robust objects the posterior median solution is quiescent, but we cannot confidently exclude star-forming solutions. Additional columns for this table, including spectral fitting results from Eazy and LePhare (see Section 3.3) are provided as supplementary online material. ID | RA | DEC | F150W | F200W | Redshift |  $t_{\mathrm{form}}$ / Gyr | $z_{\mathrm{form}}$ | log${}_{10}(M_{*}/$M${}_{\odot})$ | $U-V$ | $V-J$ | Robust
---|---|---|---|---|---|---|---|---|---|---|---
17318 | 214.808153 | 52.832201 | 27.33 | 25.73 | $4.50^{+0.13}_{-0.10}$ | $0.5^{+0.3}_{-0.2}$ | $9.5^{+6.2}_{-2.9}$ | $10.13^{+0.04}_{-0.04}$ | $1.22^{+0.13}_{-0.09}$ | $0.26^{+0.14}_{-0.15}$ | True
28316 | 214.871228 | 52.845073 | 23.44 | 22.22 | $3.53^{+0.12}_{-0.12}$ | $1.6^{+0.1}_{-0.1}$ | $3.8^{+0.2}_{-0.1}$ | $10.84^{+0.04}_{-0.04}$ | $1.06^{+0.03}_{-0.04}$ | $0.57^{+0.12}_{-0.12}$ | False
29497 | 214.760622 | 52.845322 | 22.95 | 21.57 | $3.25^{+0.08}_{-0.08}$ | $1.2^{+0.2}_{-0.3}$ | $4.9^{+1.4}_{-0.6}$ | $11.34^{+0.06}_{-0.06}$ | $1.44^{+0.10}_{-0.06}$ | $0.53^{+0.10}_{-0.08}$ | True
36262 | 214.895614 | 52.856497 | 22.77 | 21.68 | $3.26^{+0.09}_{-0.10}$ | $1.2^{+0.2}_{-0.3}$ | $4.8^{+1.2}_{-0.5}$ | $11.06^{+0.12}_{-0.07}$ | $1.11^{+0.10}_{-0.04}$ | $0.11^{+0.09}_{-0.08}$ | True
40015 | 214.853899 | 52.861358 | 25.25 | 23.23 | $3.68^{+0.16}_{-0.20}$ | $1.1^{+0.2}_{-0.2}$ | $5.2^{+1.1}_{-0.6}$ | $11.53^{+0.10}_{-0.10}$ | $2.04^{+0.08}_{-0.08}$ | $1.44^{+0.16}_{-0.15}$ | False
42128 | 214.850568 | 52.866030 | 25.35 | 23.81 | $4.19^{+0.12}_{-0.12}$ | $0.4^{+0.3}_{-0.2}$ | $10.8^{+8.2}_{-3.4}$ | $11.06^{+0.07}_{-0.05}$ | $1.48^{+0.14}_{-0.07}$ | $0.65^{+0.13}_{-0.12}$ | True
44362 | 214.879163 | 52.869187 | 25.32 | 24.48 | $3.39^{+0.29}_{-0.14}$ | $1.7^{+0.1}_{-0.2}$ | $3.7^{+0.3}_{-0.1}$ | $9.63^{+0.04}_{-0.03}$ | $0.78^{+0.03}_{-0.03}$ | $0.07^{+0.09}_{-0.08}$ | False
52124 | 214.866027 | 52.884091 | 24.7 | 23.23 | $3.38^{+0.19}_{-0.09}$ | $1.2^{+0.3}_{-0.3}$ | $4.8^{+1.2}_{-0.8}$ | $10.81^{+0.11}_{-0.05}$ | $1.52^{+0.05}_{-0.09}$ | $0.76^{+0.14}_{-0.17}$ | False
52175 | 214.866039 | 52.884255 | 23.84 | 22.31 | $3.44^{+0.14}_{-0.08}$ | $1.2^{+0.1}_{-0.2}$ | $4.8^{+0.7}_{-0.3}$ | $10.87^{+0.05}_{-0.02}$ | $1.32^{+0.04}_{-0.10}$ | $0.33^{+0.05}_{-0.10}$ | True
75768 | 214.904841 | 52.935352 | 24.51 | 23.25 | $3.31^{+0.10}_{-0.07}$ | $1.3^{+0.2}_{-0.3}$ | $4.6^{+1.2}_{-0.4}$ | $10.48^{+0.07}_{-0.04}$ | $1.25^{+0.07}_{-0.06}$ | $0.30^{+0.08}_{-0.09}$ | True
80785 | 214.915559 | 52.949026 | 24.76 | 23.96 | $4.86^{+0.15}_{-0.16}$ | $0.5^{+0.2}_{-0.2}$ | $9.2^{+3.6}_{-2.1}$ | $10.90^{+0.08}_{-0.05}$ | $1.04^{+0.09}_{-0.04}$ | $0.08^{+0.18}_{-0.08}$ | False
8888 | 214.767258 | 52.817698 | 25.19 | 23.73 | $3.49^{+0.17}_{-0.12}$ | $1.2^{+0.2}_{-0.3}$ | $4.9^{+1.1}_{-0.6}$ | $10.54^{+0.12}_{-0.08}$ | $1.34^{+0.06}_{-0.06}$ | $0.43^{+0.12}_{-0.10}$ | True
92564 | 214.957874 | 52.980293 | 24.9 | 23.35 | $3.47^{+0.11}_{-0.10}$ | $1.2^{+0.1}_{-0.1}$ | $4.9^{+0.5}_{-0.4}$ | $10.49^{+0.08}_{-0.04}$ | $1.29^{+0.07}_{-0.07}$ | $0.30^{+0.09}_{-0.11}$ | True
97581 | 214.981800 | 52.991238 | 24.14 | 22.64 | $3.46^{+0.10}_{-0.09}$ | $1.0^{+0.2}_{-0.2}$ | $5.5^{+1.0}_{-0.8}$ | $10.81^{+0.06}_{-0.03}$ | $1.35^{+0.07}_{-0.10}$ | $0.40^{+0.05}_{-0.11}$ | True
101962 | 215.039054 | 53.002778 | 26.39 | 24.85 | $4.39^{+0.11}_{-0.14}$ | $0.4^{+0.3}_{-0.2}$ | $12.1^{+6.3}_{-4.5}$ | $10.63^{+0.06}_{-0.04}$ | $1.36^{+0.11}_{-0.07}$ | $0.46^{+0.16}_{-0.11}$ | True

###  3.3 Comparison of our results with other codes

As an additional check on our results, we fit our photometry for these 15
objects with two additional SED fitting codes. We firstly run Eazy (Brammer et
al., 2008), using both the Pegase and principal component analysis (PCA)
template sets, and additionally run LePhare (Arnouts et al., 1999; Ilbert et
al., 2006), using the Bruzual & Charlot (2003) models. The redshifts for these
three spectral fitting runs, along with best-fit masses and sSFRs from LePhare
are provided in the online version of Table 2 as supplementary material.

For our robust sub-sample, all three additional sets of redshifts are in good
agreement with our Bagpipes results, with typical variations of d$z\simeq
0.1-0.2$. The masses returned by LePhare are also very similar to our Bagpipes
posterior median values, with a mean offset of 0.07 dex. LePhare also returns
sSFRs below the threshold defined in Equation 1 for all 10 objects in our
robust sub-sample.

For our 5 non-robust objects, the agreement is similar for 4 objects, however
for the remaining object (ID: 44362), LePhare finds a lower-mass star-forming
solution, though all of the codes still find this object to be at $z\simeq 3$.
This finding reflects the lower level of certainty we attach to objects that
do not meet our robust selection criteria.

##  4 Results

From our analysis, we identify a total of 15 quiescent galaxies at $z>3$, of
which 10 are members of our robust sub-sample. From these 15 objects, 4 are
placed reliably at $z>4$, and 3 of these are robustly identified as quiescent.
This is the first clear identification of massive quiescent galaxies
significantly beyond $z=4$. We provide coordinates, photometric redshifts,
magnitudes and physical properties for the 15 quiescent galaxies we identify
in Table 2. We present SEDs for our 10 robust objects in Figs 1 and 7. Cutout
images for these 10 robust galaxies are presented in Fig. 8.

###  4.1 Robust massive quiescent galaxies at $\mathbf{z>4}$

SEDs and colour images for the three $z>4$ objects in our robust sub-sample
are shown in Fig. 1. All three of these display a strong Balmer break between
the F200W and F277W bands, which provides a very strong constraint on their
redshifts. They also exhibit well-constrained red spectral slopes in the rest-
frame near-UV, indicating a lack of ongoing star formation, and blue spectral
slopes in the longer-wavelength NIRCam bands, strongly ruling our lower-
redshift dusty solutions. In this section we briefly discuss the observed
properties of each object, before moving on to discuss their SFHs in Section
5.

####  4.1.1 Galaxy 17318

The highest-redshift robust quiescent galaxy in our sample is object 17318,
with F200W = 25.7 and a photometric redshift of $z~{}\simeq~{}4.5$. This
galaxy, as shown in the top panel of Fig. 1, exhibits the characteristic
triangular SED shape of post-starburst (PSB) galaxies at lower redshift (e.g.
Wild et al. 2014), implying a recent, rapid fall in SFR (e.g. Wild et al.
2020; D’Eugenio et al. 2021). This is perhaps unsurprising, considering this
object is observed just $\simeq 1.4$ billion years after the Big Bang. For
this galaxy we derive a stellar mass of
log${}_{10}(M_{*}/$M${}_{\odot})=10.13\pm 0.06$, and a $2\sigma$ upper
limiting sSFR (97.5th percentile) of log${}_{10}($sSFR/yr${}^{-1})=-10.6$.

Our sample is shown on the rest-frame UVJ colour diagram in Fig. 3. Galaxy
17318 is the bluest (closest to the bottom-left) robust object shown in the
central panel, just below the horizontal edge of the solid UVJ selection box,
again highly consistent with lower-redshift PSBs (e.g. Belli et al. 2019;
Carnall et al. 2019a). The locations of $z>3$ quiescent galaxies on various
other rest-frame colour selection diagrams (e.g., NUVrJ; Ilbert et al. 2013)
will be explored in upcoming work by Gould et al. in prep.

This object is at a similar redshift to GOODSS-9209, the highest-redshift
candidate identified by Carnall et al. (2020), which is the target of NIRSpec
Cycle 1 observations (Carnall et al., 2023a). This new galaxy however is
$\simeq 0.6$ dex less massive, and approximately 2 magnitudes fainter at
$\lambda=2\,\mu$m (GOODSS-9209 has $K_{s}=23.6$).

This object is not included in the CANDELS catalogue of Stefanon et al.
(2017). We measure a F150W magnitude for this object of 27.4, which is fainter
than the 50 per cent point-source completeness limit calculated by Stefanon et
al. (2017) in their F160W selection band of 27.23. It is therefore not
surprising that this object was not included.

####  4.1.2 Galaxy 101962

The second-highest-redshift robust candidate we identify is object 101692 at
$z~{}\simeq~{}4.4$, shown in the middle panel of Fig. 1. It is slightly redder
than 17318, and is the middle robust quiescent object shown on the UVJ diagram
in the central panel of Fig. 3. In most respects however, this object is quite
similar to our other $z>4$ robust candidates. For this galaxy we derive a
stellar mass of log${}_{10}(M_{*}/$M${}_{\odot})=10.63^{+0.06}_{-0.04}$
(approximately 0.5 dex more massive than 17318) and a $2\sigma$ upper limiting
sSFR of log${}_{10}($sSFR/yr${}^{-1})=-9.9$.

This is the only one of our three robust $z>4$ candidates to appear in the
CANDELS EGS catalogue of Stefanon et al. (2017), having a match within 0.1′′
(ID: 23297). We show the CANDELS photometry for this object in blue in Fig. 2.
This galaxy was only previously detected with $3\sigma$ significance in the
WFCAM $K_{s}$ band, and no data were available at all between
$\lambda=2.3-3.0\mu$m, meaning virtually no constraint could be placed on the
Balmer break. The median photometric redshift reported by the CANDELS team
(following the method of Dahlen et al. 2013) is $z=3.5$. However, the
individual estimates from their different codes display considerable variance,
with a $1\sigma$ range from $z=2.45-4.70$.

As a demonstration, we fit the Stefanon et al. (2017) photometry shown in Fig.
2 with Bagpipes, using the same model described in Section 3.1. Our best fit
is shown in green in Fig. 2, along with the gold JWST data and red best fit
from Fig. 1. From the Stefanon et al. (2017) data we recover an extremely
broad photometric redshift posterior, $z=1.34^{+1.94}_{-0.39}$. The 1D and 2D
stellar mass, sSFR and redshift posteriors for both fits are shown in the
corner plot on the left side of Fig. 2. It can be seen that the lack of a
well-constrained redshift results in a very broad stellar mass posterior, and
virtually no constraint on the sSFR, with the sSFR posterior following the
prior imposed by our SFH model (e.g., see Fig. 1 of Carnall et al. 2019b).

####  4.1.3 Galaxy 42128

Galaxy 42128 has a slightly lower redshift, $z\simeq 4.2$, and is both
brighter and more massive than the other two galaxies. It is also the reddest
of our $z>4$ candidates on the UVJ diagram. However, in most respects, this
galaxy is very similar to 17318 and 101962, also exhibiting the characteristic
triangular PSB spectral shape. For this galaxy, we derive a stellar mass of
log${}_{10}(M_{*}/$M${}_{\odot})=11.06^{+0.07}_{-0.05}$, and a $2\sigma$ upper
limiting sSFR of log${}_{10}($sSFR/yr${}^{-1})=-11.2$.

In shorter-wavelength, lower spatial resolution imaging this object is not
distinguishable from the extended structure of the nearby barred-spiral galaxy
and therefore does not feature in the Stefanon et al. (2017) CANDELS EGS
catalogue. However, the longer wavelength and higher spatial resolution
imaging provided by JWST reveal this object as a brighter, redder, highly
compact background source.

Table 3: Number densities derived for our quiescent galaxy sample in integer
redshift bins spanning $3<z<5$. Uncertainties were calculated as the Poisson
noise on the number of objects found (Gehrels, 1986). The full quiescent
sample includes quiescent galaxies for which we cannot rule out secondary
star-forming solutions, whereas the robust sub-sample includes only those
galaxies for which we can confidently exclude star-forming solutions. Full
quiescent sample
---
| $N_{\mathrm{galaxies}}$ |  $n$ / Mpc-3 |  $n$ / Mpc-3
Redshift range | F200W < 26.5 | F200W < 26.5 | F200W < 24.5
$3<z<4$ | 11 | $11.6^{+4.7}_{-3.4}\times 10^{-5}$ | $10.6^{+4.5}_{-3.3}\times 10^{-5}$
$4<z<5$ | 4 | $4.7^{+3.7}_{-2.2}\times 10^{-5}$ | $2.3^{+3.1}_{-1.5}\times 10^{-5}$
Robust sub-sample
Redshift range | $N_{\mathrm{galaxies}}$ |  $n$ / Mpc-3 |  $n$ / Mpc-3
| F200W < 26.5 | F200W < 26.5 | F200W < 24.5
$3<z<4$ | 7 | $7.4^{+4.0}_{-2.7}\times 10^{-5}$ | $6.3^{+3.8}_{-2.5}\times 10^{-5}$
$4<z<5$ | 3 | $3.5^{+3.4}_{-1.9}\times 10^{-5}$ | $1.2^{+2.7}_{-1.0}\times 10^{-5}$
Figure 4: Number density estimates for high-redshift massive quiescent
galaxies. Our estimate at $3<z<5$ derived from the JWST CEERS data are $3-5$
times higher than pre-JWST estimates, and, at $z\simeq 3$, approach the result
of McLeod et al. (2021) for the total galaxy population. Stellar masses
derived by other authors have been converted to a Kroupa (2001) IMF where
necessary.

###  4.2 Comparison with previous work in the EGS field

A previous search for $z>3$ massive quiescent galaxies in the EGS field was
presented by Merlin et al. (2019), who analysed the CANDELS catalogue of
Stefanon et al. (2017). In this subsection we present a comparison with their
results to demonstrate the power of the new JWST data. From our 15-object
sample, 12 objects have matches in the Stefanon et al. (2017) catalogue,
including object 101963 discussed in Section 4.1.2.

Merlin et al. (2019) identify 13 candidate $z>3$ quiescent galaxies in the
Stefanon et al. (2017) catalogue. From these, 4 objects fall within the area
covered by our JWST CEERS catalogue (a further object falls within a gap
between two NIRCam short-wavelength detectors). All 4 of these are also
included in our sample (IDs: 28316, 36262, 52175, 75768). These are all at
$z\lesssim 3.5$, and are the 4 brightest galaxies in our sample in the F150W
band. The other 8 objects in our sample that also appear in the Stefanon et
al. (2017) catalogue were not selected by Merlin et al. (2019), likely as a
result of their relative faintness, which makes constraining redshifts, masses
and sSFRs more challenging, as demonstrated in Fig. 2.

###  4.3 Quiescent galaxy number densities at $\mathbf{3<z<5}$

In Table 3 we report our estimates of the number densities of quiescent
galaxies over the redshift range $3<z<5$ that meet our F200W < 26.5 selection
threshold. We report number densities based on both our full quiescent sample
and our robust sub-sample. The numbers based on the robust sub-sample can be
interpreted as conservative lower limits.

In Schreiber et al. (2018), the authors report a number density of $2.0\pm
0.3\ \times\ 10^{-5}$ Mpc-3 for a spectroscopic sample of quiescent galaxies
at $3<z<4$ with $K_{s}<24.5$. We compare our results with Schreiber et al.
(2018) by calculating the number of F200W < 24.5 quiescent galaxies in our
sample at $3<z<4$, as this filter is closest in wavelength coverage to the
$K_{s}$ band. These numbers are reported in the right-hand column of Table 3.

The number density we find for our robust sub-sample at F200W < 24.5 is
approximately a factor of 3 larger than the result of Schreiber et al. (2018),
whereas the result from our full quiescent sample is approximately a factor of
5 larger. We attribute this significant increase in the number of $z>3$
quiescent galaxies to the much deeper, redder imaging now available from JWST,
which allows physical properties to be reliably inferred for faint, red
galaxies such as these.

From our 11 quiescent galaxies at $3<z<4$, a total of 10 have matches within
0.25′′ in the Stefanon et al. (2017) CANDELS EGS catalogue. From these 10, a
total of 9 have CANDELS photometric redshifts in the range from $3<z<4$, and
by fitting the Stefanon et al. (2017) photometry with Bagpipes we recover
similar results. However, just 3 of these objects would be included in our
quiescent sample given the CANDELS photometry, as their sSFRs are far less
well constrained by these data. Just 2 of these 3 would be identified as
robust.

It is therefore likely that, had we designed a spectroscopic follow-up
campaign similar to that of Schreiber et al. (2018) based on fitting only
CANDELS EGS data, we would have arrived at a very similar number density for
$3<z<4$ quiescent galaxies to the one they obtain. Indeed, our result based on
fitting CANDELS UDS and GOODS South photometry in Carnall et al. (2020)
arrived at a number density of $1.7\pm 0.3\ \times\ 10^{-5}$ Mpc-3 for $3<z<4$
quiescent galaxies with $K_{s}<24.5$, fully consistent with Schreiber et al.
(2018).

In Fig. 4 we show a comparison of $z>3$ quiescent galaxy number density
estimates from the literature, including our new results. We restrict this
comparison to galaxies with log${}_{10}(M_{*}/\mathrm{M_{\odot}})>10.5$, where
the CANDELS catalogues are mass complete at $z<4$. At $z>4$, the Balmer break
moves to $\lambda>2\mu$m, meaning selection at longer wavelengths is likely to
be necessary to obtain mass-complete samples.

As above, it can be seen from Fig. 4 that our new result (for our full sample)
is a factor of $\simeq 5$ higher than previous estimates, which are all in
close agreement. We also show results at $z<3.75$ derived from the stellar
mass functions of McLeod et al. (2021), both for quiescent galaxies and for
the whole galaxy population. Our $3<z<4$ quiescent galaxy number density is
close to the total number density found by McLeod et al. (2021). This suggests
the small area studied in this work may be over-dense at these redshifts.

Finally we note that, as shown in Carnall et al. (2019b), standard parametric
SFH models impose strong priors on galaxy sSFRs, favouring sSFRs close to the
star-forming main sequence. In the absence of strongly constraining data, this
could plausibly result in quiescent galaxies being misidentified as star
forming. This effect is demonstrated in Fig. 2, in which the sSFR posterior
for galaxy 101962 using CANDELS data is strongly weighted towards star-forming
solutions, even at the correct redshift of $z\simeq 4.4$. Clearly, follow-up
studies using larger-area JWST imaging surveys are of critical importance to
clarify this situation.

Figure 5: Formation redshifts for $z>3$ massive quiescent galaxy sample, shown
as a function of their observed redshifts. Our candidates at $z>4$ have
formation redshifts from $9<z_{\mathrm{form}}<12$, whereas our $3<z<4$
galaxies are all younger, having formed at $z_{\mathrm{form}}<6$. Figure 6:
Star-formation histories for our 10 robustly identified $z>3$ quiescent
galaxies. The solid orange lines show the 50th percentiles of the SFH
posterior distributions in each case, whereas the orange shaded regions show
the $1\sigma$ range. The dashed gray vertical lines show the posterior median
redshift of formation for each object. SEDs are shown in Figs 1 and 7.

##  5 Discussion

In this section we discuss the star-formation histories we infer for our
sample. We calculate formation times, $t_{\mathrm{form}}$, for each object as
the average time at which the stars in the galaxy formed. This is the age of
the Universe at the time corresponding to the (mass-weighted) mean stellar age
(see Carnall et al. 2018, Equation 11). We also calculate formation redshifts,
$z_{\mathrm{form}}$, which are the redshifts corresponding to
$t_{\mathrm{form}}$. These values are reported in Table 2.

Formation redshifts for our sample are plotted against the observed redshift
of each galaxy in Fig. 5. In addition, the full SFH posteriors for the 10
galaxies in our robust sub-sample are shown in Fig. 6. Their posterior median
formation redshifts are shown with dashed gray vertical lines.

It can be seen that the 3 robust quiescent galaxies at $z>4$ formed the bulk
of their stellar populations during the first billion years at $z>6$, with
posterior median formation redshifts in the range $9<z_{\mathrm{form}}<12$.
The second-highest-redshift object, 101962, is the earliest formed, with the
bulk of its stars having formed at $z>8$. These findings make our $z>4$
objects highly plausible as descendants of the sample identified by Labbe et
al. (2022). Indeed, objects 101962 and 42128 are both predicted to have had
log${}_{10}(M_{*}/$M${}_{\odot})>10$ at their formation redshifts of
$z_{\mathrm{form}}=12.1^{+6.3}_{-4.5}$ and
$z_{\mathrm{form}}=10.8^{+8.2}_{-3.4}$ respectively. This is fully consistent
with the finding of Labbe et al. (2022) that a considerable number of
log${}_{10}(M_{*}/$M${}_{\odot})>10$ galaxies were already in place by
$7<z<11$.

The 7 robust quiescent galaxies at $3<z<4$ formed their stellar populations
later in cosmic time, with $4<z_{\mathrm{form}}<6$. The fact that none of
these galaxies is older suggests that the $z>4$ quiescent galaxies in our
sample have not yet reached the end-point of their evolution, and are likely
to experience further star formation by $z\simeq 3$. However, larger-area JWST
surveys and spectroscopic follow-up will certainly be required to rule out the
possibility of $3<z<4$ quiescent galaxies with stellar populations dating back
to $z>6$.

It is interesting to reflect on the fact that we do not identify any quiescent
galaxies at $z>5$, despite the fact that our oldest object appears to have
quenched by $z>6$. Apart from the obvious limitations of the relatively small
imaging area included in the initial CEERS release, it is currently unclear
how the colours of newly quenched galaxies evolve in order to arrive in the
UVJ-quiescent box, and how long this might take after star formation ceases
(e.g. Belli et al. 2019; Carnall et al. 2019b; Akins et al. 2022).

Upcoming larger-area JWST imaging surveys, such as Public Release Imaging for
Extragalactic Research (PRIMER333https://primer-jwst.github.io) are ideally
suited to searching for evidence of quiescent galaxies at $z>5$, as well as
selecting larger samples of $3<z<5$ quiescent galaxies to produce more-robust
number densities. Such searches may benefit from the use of catalogues
selected in longer-wavelength NIRCam bands, rather than the F200W selection
employed in this work, as well as new, JWST-specific colour selection criteria
(e.g. Leja et al. 2019; Antwi-Danso et al. 2022).

##  6 Conclusion

In this work we present the results of a search for massive quiescent galaxies
at redshifts $z>3$, selected from the first NIRCam data taken by the JWST
CEERS Early Release Science programme. We identify 15 galaxies in the redshift
range $3<z<5$ with robust photometric redshifts and posterior median sSFRs
that suggest they are quiescent. For 10 of these galaxies, we can confidently
rule out star-forming solutions, and we refer to these galaxies as comprising
our robust sub-sample.

Three of our robustly quiescent objects are at $z>4$, and these constitute the
best evidence to date for the existence of massive quiescent galaxies
significantly above $z=4$. These objects have stellar masses in the range
log${}_{10}(M_{*}/$M${}_{\odot})=10.1-11.1$ and posterior median formation
redshifts from $9<z_{\mathrm{form}}<12$. Two of these $z>4$ galaxies would
have had stellar masses in excess of log${}_{10}(M_{*}/$M${}_{\odot})=10$ by
$z\gtrsim 8$, supporting the recent findings of Labbe et al. (2022).

Only one of these $z>4$ objects has a match in the Stefanon et al. (2017)
CANDELS EGS catalogue. The redshift of this object was not well constrained by
previously available data, with both the CANDELS team and Bagpipes returning
extremely broad redshift posteriors, with best-fitting values at $z<4$.

We calculate number densities for our quiescent sample, as well as for our
robust sub-sample (which can be regarded as yielding robust lower limits). We
find that the number density of $3<z<4$ quiescent galaxies with F200W < 24.5
is a factor of $3-5$ times higher than previously reported by Schreiber et al.
(2018). We demonstrate that this difference arises as a result of better
constraints from the new NIRCam data, and show that previously available data
from CANDELS would have led us to a similar number density to that calculated
by Schreiber et al. (2018). This finding poses an additional challenge for
simulations of early galaxy formation, which already struggle to reproduce
previously reported number densities.

The 10 robust $z>3$ massive quiescent galaxies we report, and the 3 at $z>4$
in particular, are excellent potential targets for follow up NIRSpec
spectroscopy, either as part of CEERS, or via dedicated programmes in JWST
Cycle 2 and beyond. The detection of Balmer absorption features would
unambiguously confirm these objects as quiescent at $z>4$, and full spectral
fitting of deep continuum spectroscopic data would provide strong constraints
on their SFHs, as well as the ability to probe in detail their physical
properties.

## Acknowledgements

A. C. Carnall thanks the Leverhulme Trust for their support via a Leverhulme
Early Career Fellowship. R. Begley, D. J. McLeod, M. L. Hamadouche, C. Donnan,
R. J. McLure, J. S. Dunlop and V. Wild acknowledge the support of the Science
and Technology Facilities Council. F. Cullen acknowledges support from a UKRI
Frontier Research Guarantee Grant [grant reference EP/X021025/1]. S. Jewell
and C. Pollock acknowledge the support of the School of Physics & Astronomy,
University of Edinburgh via Summer Studentship bursaries.

## Data Availability

All JWST and HST data products are available via the Mikulski Archive for
Space Telescopes (https://mast.stsci.edu). Photometric data and fitted model
posteriors are available upon request.

## References

  * Akins et al. (2022) Akins H. B., Narayanan D., Whitaker K. E., Davé R., Lower S., Bezanson R., Feldmann R., Kriek M., 2022, ApJ, 929, 94
  * Antwi-Danso et al. (2022) Antwi-Danso J., et al., 2022, arXiv e-prints, p. arXiv:2207.07170
  * Arnouts et al. (1999) Arnouts S., Cristiani S., Moscardini L., Matarrese S., Lucchin F., Fontana A., Giallongo E., 1999, MNRAS, 310, 540
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Belli et al. (2019) Belli S., Newman A. B., Ellis R. S., 2019, ApJ, 874, 17
  * Bertin (2006) Bertin E., 2006, in Gabriel C., Arviset C., Ponz D., Enrique S., eds, Astronomical Society of the Pacific Conference Series Vol. 351, Astronomical Data Analysis Software and Systems XV. p. 112
  * Bertin (2010) Bertin E., 2010, SWarp: Resampling and Co-adding FITS Images Together, Astrophysics Source Code Library, record ascl:1010.068 (ascl:1010.068)
  * Bertin & Arnouts (1996) Bertin E., Arnouts S., 1996, A&AS, 117, 393
  * Beverage et al. (2021) Beverage A. G., Kriek M., Conroy C., Bezanson R., Franx M., van der Wel A., 2021, ApJ, 917, L1
  * Brammer et al. (2008) Brammer G. B., van Dokkum P. G., Coppi P., 2008, ApJ, 686, 1503
  * Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Buchner et al. (2014) Buchner J., et al., 2014, A&A, 564, A125
  * Byler et al. (2017) Byler N., Dalcanton J. J., Conroy C., Johnson B. D., 2017, ApJ, 840, 44
  * Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, ApJ, 533, 682
  * Carnall et al. (2018) Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379
  * Carnall et al. (2019a) Carnall A. C., et al., 2019a, MNRAS, 490, 417
  * Carnall et al. (2019b) Carnall A. C., Leja J., Johnson B. D., McLure R. J., Dunlop J. S., Conroy C., 2019b, ApJ, 873, 44
  * Carnall et al. (2020) Carnall A. C., et al., 2020, MNRAS, 496, 695
  * Carnall et al. (2022) Carnall A. C., et al., 2022, ApJ, 929, 131
  * Carnall et al. (2023a) Carnall A. C., et al., 2023a, arXiv e-prints, p. arXiv:2301.11413
  * Carnall et al. (2023b) Carnall A. C., et al., 2023b, MNRAS, 518, L45
  * Castellano et al. (2022) Castellano M., et al., 2022, ApJ, 938, L15
  * Cecchi et al. (2019) Cecchi R., Bolzonella M., Cimatti A., Girelli G., 2019, ApJ, 880, L14
  * Charlot & Fall (2000) Charlot S., Fall S. M., 2000, ApJ, 539, 718
  * Chevallard & Charlot (2016) Chevallard J., Charlot S., 2016, MNRAS, 462, 1415
  * D’Eugenio et al. (2021) D’Eugenio C., et al., 2021, A&A, 653, A32
  * Dahlen et al. (2013) Dahlen T., et al., 2013, ApJ, 775, 93
  * Davé et al. (2016) Davé R., Thompson R., Hopkins P. F., 2016, MNRAS, 462, 3265
  * Dickinson & FIDEL Team (2007) Dickinson M., FIDEL Team 2007, in American Astronomical Society Meeting Abstracts. p. 52.16
  * Donnan et al. (2023) Donnan C. T., et al., 2023, MNRAS, 518, 6011
  * Esdaile et al. (2021) Esdaile J., et al., 2021, AJ, 162, 225
  * Falcón-Barroso et al. (2011) Falcón-Barroso J., Sánchez-Blázquez P., Vazdekis A., Ricciardelli E., Cardiel N., Cenarro A. J., Gorgas J., Peletier R. F., 2011, A&A, 532, A95
  * Feroz et al. (2019) Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2019, The Open Journal of Astrophysics, 2, 10
  * Finkelstein et al. (2022) Finkelstein S. L., et al., 2022, ApJ, 940, L55
  * Forrest et al. (2020) Forrest B., et al., 2020, ApJ, 903, 47
  * Gallazzi et al. (2014) Gallazzi A., Bell E. F., Zibetti S., Brinchmann J., Kelson D. D., 2014, ApJ, 788, 72
  * Gehrels (1986) Gehrels N., 1986, ApJ, 303, 336
  * Girelli et al. (2019) Girelli G., Bolzonella M., Cimatti A., 2019, A&A, 632, A80
  * Glazebrook et al. (2017) Glazebrook K., et al., 2017, Nature, 544, 71
  * Hamadouche et al. (2022) Hamadouche M. L., et al., 2022, MNRAS, 512, 1262
  * Ilbert et al. (2006) Ilbert O., et al., 2006, A&A, 457, 841
  * Ilbert et al. (2013) Ilbert O., et al., 2013, A&A, 556, A55
  * Inoue et al. (2014) Inoue A. K., Shimizu I., Iwata I., Tanaka M., 2014, MNRAS, 442, 1805
  * Koekemoer et al. (2011) Koekemoer A. M., et al., 2011, The Astrophysical Journal Supplement Series, 197, 36
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Labbe et al. (2022) Labbe I., et al., 2022, arXiv e-prints, p. arXiv:2207.12446
  * Leja et al. (2019) Leja J., Tacchella S., Conroy C., 2019, ApJ, 880, L9
  * Marigo et al. (2013) Marigo P., Bressan A., Nanni A., Girardi L., Pumo M. L., 2013, MNRAS, 434, 488
  * Marsan et al. (2022) Marsan Z. C., et al., 2022, ApJ, 924, 25
  * McLeod et al. (2016) McLeod D. J., McLure R. J., Dunlop J. S., 2016, MNRAS, 459, 3812
  * McLeod et al. (2021) McLeod D. J., McLure R. J., Dunlop J. S., Cullen F., Carnall A. C., Duncan K., 2021, MNRAS, 503, 4413
  * Merlin et al. (2018) Merlin E., et al., 2018, MNRAS, 473, 2098
  * Merlin et al. (2019) Merlin E., et al., 2019, MNRAS, 490, 3309
  * Naidu et al. (2022) Naidu R. P., et al., 2022, ApJ, 940, L14
  * Pacifici et al. (2016) Pacifici C., et al., 2016, ApJ, 832, 79
  * Rigby et al. (2022) Rigby J., et al., 2022, arXiv e-prints, p. arXiv:2207.05632
  * Salim et al. (2018) Salim S., Boquien M., Lee J. C., 2018, ApJ, 859, 11
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez P., et al., 2006, MNRAS, 371, 703
  * Santini et al. (2021) Santini P., et al., 2021, A&A, 652, A30
  * Schreiber et al. (2018) Schreiber C., et al., 2018, A&A, 618, A85
  * Skilling (2006) Skilling J., 2006, Bayesian Anal., 1, 833
  * Stefanon et al. (2017) Stefanon M., et al., 2017, ApJS, 229, 32
  * Stevans et al. (2021) Stevans M. L., et al., 2021, ApJ, 921, 58
  * Suess et al. (2022) Suess K. A., et al., 2022, ApJ, 937, L33
  * Valentino et al. (2020) Valentino F., et al., 2020, ApJ, 889, 93
  * Wild et al. (2014) Wild V., et al., 2014, MNRAS, 440, 1880
  * Wild et al. (2020) Wild V., et al., 2020, MNRAS, 494, 529
  * Williams et al. (2009) Williams R. J., Quadri R. F., Franx M., van Dokkum P., Labbé I., 2009, ApJ, 691, 1879
  * Wu et al. (2018) Wu P.-F., et al., 2018, ApJ, 868, 37
  * van der Wel et al. (2014) van der Wel A., et al., 2014, ApJ, 788, 28

##  Appendix A SEDs and cutout images for robust quiescent galaxies

In Fig. 7 we show SEDs for the 7 robust quiescent galaxies we identify at
$3<z<4$. SEDs for the 3 robust quiescent galaxies we identify at $z>4$ are
shown in Fig. 1. In Fig. 8 we show $5^{\prime\prime}\ \times\
5^{\prime\prime}$ HST ACS + JWST NIRCam cutout images for each of the 10
galaxies in our robust sample.

Figure 7: Spectral energy distributions for our 7 robust $3<z<4$ quiescent
galaxies (SEDs for the three robust galaxies at $z>4$ are shown in Fig. 1).
Our 10-band photometric data from HST ACS and JWST NIRCam are shown in blue.
The posterior median Bagpipes models are overlaid in red. The inset panels
show the position of each object on the UVJ diagram. Figure 8: Cutout images
for each of the 10 objects in our robust quiescent sample, in descending order
of redshift. Each cutout image is $5^{\prime\prime}\ \times\
5^{\prime\prime}$.