# JADES: Rest-frame UV-to-NIR Size Evolution of Massive Quiescent Galaxies
from Redshift $z=5$ to $z=0.5$

Zhiyuan Ji Steward Observatory, University of Arizona, 933 N. Cherry Avenue,
Tucson, AZ 85721, USA  Christina C. Williams NSF’s National Optical-Infrared
Astronomy Research Laboratory, 950 N. Cherry Avenue, Tucson, AZ 85719, USA
Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ
85721, USA  Katherine A. Suess Department of Astronomy and Astrophysics,
University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064,
USA  Kavli Institute for Particle Astrophysics and Cosmology and Department of
Physics, Stanford University, Stanford, CA 94305, USA  Sandro Tacchella Kavli
Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK  Cavendish Laboratory, University of Cambridge, 19 JJ Thomson
Avenue, Cambridge, CB3 0HE, UK  Benjamin D. Johnson Center for Astrophysics
$|$ Harvard & Smithsonian, 60 Garden St., Cambridge MA 02138 USA  Brant
Robertson Department of Astronomy and Astrophysics, University of California,
Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA  Stacey Alberts
Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ
85721, USA  William M. Baker Kavli Institute for Cosmology, University of
Cambridge, Madingley Road, Cambridge, CB3 0HA, UK  Cavendish Laboratory,
University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK  Stefi
Baum 2Department of Physics and Astronomy, University of Manitoba, Winnipeg,
MB R3T 2N2, Canada  Rachana Bhatawdekar European Space Agency (ESA), European
Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, 28692 Villanueva
de la Cañada, Madrid, Spain  Nina Bonaventura Steward Observatory, University
of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA  Kristan Boyett School
of Physics, University of Melbourne, Parkville 3010, VIC, Australia  ARC
Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D),
Australia  Andrew J. Bunker Department of Physics, University of Oxford, Denys
Wilkinson Building, Keble Road, Oxford OX1 3RH, UK  Stefano Carniani Scuola
Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy  Stephane
Charlot Sorbonne Université, CNRS, UMR 7095, Institut d’Astrophysique de
Paris, 98 bis bd Arago, 75014 Paris, France  Zuyi Chen Steward Observatory,
University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA  Jacopo
Chevallard Department of Physics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford OX1 3RH, UK  Emma Curtis-Lake Centre for
Astrophysics Research, Department of Physics, Astronomy and Mathematics,
University of Hertfordshire, Hatfield AL10 9AB, UK  Francesco D’Eugenio Kavli
Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK  Cavendish Laboratory, University of Cambridge, 19 JJ Thomson
Avenue, Cambridge, CB3 0HE, UK  Anna de Graaff Max-Planck-Institut für
Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany  Christa DeCoursey
Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ
85721, USA  Eiichi Egami Steward Observatory, University of Arizona, 933 N.
Cherry Avenue, Tucson, AZ 85721, USA  Daniel J. Eisenstein Center for
Astrophysics $|$ Harvard & Smithsonian, 60 Garden St., Cambridge MA 02138 USA
Kevin Hainline Steward Observatory, University of Arizona, 933 N. Cherry
Avenue, Tucson, AZ 85721, USA  Ryan Hausen Department of Physics and
Astronomy, The Johns Hopkins University, 3400 N. Charles St., Baltimore, MD
21218  Jakob M. Helton Steward Observatory, University of Arizona, 933 N.
Cherry Avenue, Tucson, AZ 85721, USA  Tobias J. Looser Kavli Institute for
Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue,
Cambridge, CB3 0HE, UK  Jianwei Lyu Steward Observatory, University of
Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA  Roberto Maiolino Kavli
Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK  Cavendish Laboratory, University of Cambridge, 19 JJ Thomson
Avenue, Cambridge, CB3 0HE, UK  Department of Physics and Astronomy,
University College London, Gower Street, London WC1E 6BT, UK  Michael V.
Maseda Department of Astronomy, University of Wisconsin-Madison, 475 N.
Charter St., Madison, WI 53706 USA  Erica Nelson Department for Astrophysical
and Planetary Science, University of Colorado, Boulder, CO 80309, USA  George
Rieke Steward Observatory and Dept of Planetary Sciences, University of
Arizona 933 N. Cherry Avenue Tucson AZ 85721, USA  Marcia Rieke Steward
Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721,
USA  Hans-Walter Rix Max-Planck-Institut f”̈ur Astronomie, K”̈onigstuhl 17,
D-69117, Heidelberg, Germany  Lester Sandles Kavli Institute for Cosmology,
University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK  Cavendish
Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE,
UK  Fengwu Sun Steward Observatory, University of Arizona, 933 N. Cherry
Avenue, Tucson, AZ 85721, USA  Hannah Übler Kavli Institute for Cosmology,
University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK  Cavendish
Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE,
UK  Christopher N. A. Willmer Steward Observatory, University of Arizona, 933
N. Cherry Avenue, Tucson, AZ 85721, USA  Chris Willott NRC Herzberg, 5071 West
Saanich Rd, Victoria, BC V9E 2E7, Canada  Joris Witstok Kavli Institute for
Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue,
Cambridge, CB3 0HE, UK

###### Abstract

We present the UV-to-NIR size evolution of a sample of 161 quiescent galaxies
with $M_{*}>10^{10}M_{\sun}$ over $0.5<z<5$. With deep multi-band NIRCam
images in GOODS-South from JADES, we measure the effective radii ($R_{e}$) of
the galaxies at rest-frame 0.3, 0.5 and 1$\micron$. On average, we find that
quiescent galaxies are 45% (15%) more compact at rest-frame 1$\micron$ than
they are at 0.3$\micron$ (0.5$\micron$). Regardless of wavelengths, the
$R_{e}$ of quiescent galaxies strongly evolves with redshift, and this
evolution depends on stellar mass. For lower-mass quiescent galaxies with
$M_{*}=10^{10}-10^{10.6}M_{\sun}$, the evolution follows
$R_{e}\propto(1+z)^{-1.1}$, whereas it becomes steeper, following
$R_{e}\propto(1+z)^{-1.7}$, for higher-mass quiescent galaxies with
$M_{*}>10^{10.6}M_{\sun}$. To constrain the physical mechanisms driving the
apparent size evolution, we study the relationship between $R_{e}$ and the
formation redshift ($z_{\rm{form}}$) of quiescent galaxies. For lower-mass
quiescent galaxies, this relationship is broadly consistent with
$R_{e}\propto(1+z_{\rm{form}})^{-1}$, in line with the expectation of the
progenitor effect. For higher-mass quiescent galaxies, the relationship
between $R_{e}$ and $z_{\rm{form}}$ depends on stellar age. Older quiescent
galaxies have a steeper relationship between $R_{e}$ and $z_{\rm{form}}$ than
that expected from the progenitor effect alone, suggesting that mergers and/or
post-quenching continuous gas accretion drive additional size growth in very
massive systems. We find that the $z>3$ quiescent galaxies in our sample are
very compact, with mass surface densities $\Sigma_{e}\gtrsim
10^{10}M_{\sun}/\rm{kpc}^{2}$, and their $R_{e}$ are possibly even smaller
than anticipated from the size evolution measured for lower-redshift quiescent
galaxies. Finally, we take a close look at the structure of GS-9209, one of
the earliest confirmed massive quiescent galaxies at $z_{\rm{spec}}\sim 4.7$.
From UV to NIR, GS-9209 becomes increasingly compact, and its light profile
becomes more spheroidal, showing that the color gradient is already present in
this earliest massive quiescent galaxy.

Galaxy formation(595); Galaxy evolution(594); Galaxy structure(622); High-
redshift galaxies(734)

††facilities: HST (ACS, WFC3), JWST (NIRCam)††software: Galfit (Peng et al.,
2010a), Prospector (Johnson et al., 2021), FSPS (Conroy et al., 2009; Conroy &
Gunn, 2010), MIST (Choi et al., 2016; Dotter, 2016), MILES (Falcón-Barroso et
al., 2011), astropy (Astropy Collaboration et al., 2022), photutils (Bradley
et al., 2023), Webbpsf (Perrin et al., 2014b), numpy (Harris et al., 2020),
scipy (Virtanen et al., 2020), corner (Foreman-Mackey, 2016)

##  1 Introduction

The origin and evolution of massive and quiescent galaxies that no longer form
stars are major open questions in Astrophysics. The last two decades of
research have revealed that massive quiescent galaxies emerge rapidly, but in
small numbers, at least as early as $z\gtrsim 4$. By $z<1$, such galaxies
dominate the total stellar mass budget of the Universe (e.g. Muzzin et al.,
2013; Ilbert et al., 2013). Understanding their evolutionary history is thus a
critical piece of building a coherent picture of cosmic evolution.

The near-infrared (NIR) cameras on the Hubble Space Telescope (HST) enabled
measurement of the rest-optical structure of massive galaxies up to $z\sim
2.5$. A key revelation was the extremely compact morphology of quiescent
galaxies at Cosmic Noon: They are five times smaller on average at $z\sim 1-2$
compared to their presumed descendants at $z\sim 0$ (Daddi et al., 2005;
Trujillo et al., 2007; Toft et al., 2007; van Dokkum et al., 2008; Damjanov et
al., 2011; Newman et al., 2012, just to name a few), and their stellar mass
densities rival the densest stellar structures known in the Universe (e.g.
Hopkins et al., 2010a). In addition, compared to star-forming galaxies of
similar mass and redshift, quiescent galaxies display more dramatic evolution
in size from $z\sim 2.5$ to the present (e.g., Shen et al., 2003; van der Wel
et al., 2014; Barro et al., 2017). Which processes are driving this strong
size evolution of quiescent galaxies, and whether the dominant processes
change over time have sparked contentious debate. It is expected for galaxies
to increase in size over time given the $\Lambda$CDM cosmology, because cosmic
mass density decreases over time due to the expansion of the Universe, and
galaxy size should, to some extent, reflect the density of the Universe at its
formation epoch, i.e. the progenitor effect (e.g., Peebles, 1993; Mo et al.,
2010; Lilly & Carollo, 2016). Disentangling whether the size increase towards
lower redshifts is due to the addition of larger, recently quenched sources,
or that the galaxies themselves increase in size without appreciable increase
in mass (e.g. through gas-poor minor mergers) has been a key target of
numerous studies (e.g., Bezanson et al., 2009; Hopkins et al., 2009; Naab et
al., 2009; van Dokkum et al., 2010; Fan et al., 2010; Oser et al., 2012;
Newman et al., 2012; Cassata et al., 2013; Williams et al., 2017; Suess et
al., 2019a; Ji & Giavalisco, 2022).

How the structure of quiescent galaxies evolves in earlier times – at $z>2.5$
– still eludes us. At this redshift range, HST is limited to the wavelengths
of rest-frame UV, i.e. bluewards $4000$Å. Likely because the rest-UV sizes are
not representative of the true distributions of stellar mass in galaxies due
to non-negligible impact of color gradients (e.g., Suess et al., 2019b, a;
Mosleh et al., 2020; Miller et al., 2022), earlier studies relying on HST have
found conflicting results: Some studies say that early quiescent galaxy sizes
are consistent with extrapolations from lower-redshift quiescent galaxies,
while some argue that early quiescent galaxies are exceedingly compact
(Straatman et al., 2015; Kubo et al., 2018; Lustig et al., 2021). Obviously,
in the HST era, the major limitation preventing a complete picture was the
lack of deep and high resolution rest-frame optical/NIR imaging at $z>2.5$,
which prevented both robust photometric characterization due to large
uncertainties in Spectral Energy Distribution (SED) modeling, and rest-frame
optical morphologies.

The recent launch of JWST (Gardner et al., 2023) has opened a new window to
view early quiescent galaxies in exquisite detail. Its unprecedented
sensitivity and wavelength sampling enable more robust modeling of stellar
populations and star formation history (SFH) reconstruction. The arrival of
Cycle 1 JWST data demonstrated the power of this new facility to identify and
characterize robust massive quiescent galaxy candidates at $z>2.5$, (Carnall
et al., 2023a; Valentino et al., 2023) and make spectroscopic confirmations
(Nanayakkara et al., 2022; Carnall et al., 2023b; Glazebrook et al., 2023)
even as early as $z\sim 5$. Of particular importance to understanding their
formation pathways, JWST, for the first time, now allows measurements of their
rest-frame optical and NIR structures at high angular resolution and to higher
redshifts than previously possible (out to $z\sim 10$, well beyond the limit
of HST at $z\sim 2.5$).

In this paper, we use new ultra-deep JWST/NIRCam imaging from the JWST
Advanced Deep Extragalactic Survey (JADES, Eisenstein et al. 2023) to
investigate the redshift evolution of the half-light sizes of quiescent
galaxies, measured at rest-frame UV, optical and NIR wavelengths. We focus in
particular on massive galaxies with stellar mass $\log(M_{*}/M_{\sun})>10$, in
order to focus on the formation pathways mostly independent of
environmentally-driven quenching, a process which affects primarily lower-mass
galaxies (e.g., Peng et al., 2010b; Ji et al., 2018; Papovich et al., 2018).
Throughout this paper, we adopt the AB magnitude system and the $\Lambda$CDM
cosmology with Planck Collaboration et al. 2020 parameters, i.e.,
$\Omega_{m}=0.315$ and $\rm{h=H_{0}/(100km\,s^{-1}\,Mpc^{-1})=0.673}$.

##  2 Observations

The sample of quiescent galaxies presented this work is in the GOODS-South
(Giavalisco et al., 2004) portion of JADES where we have obtained deep
JWST/NIRCam imaging in 9 filters, i.e., F090W, F115W, F150W, F200W, F277W,
F335M, F356W, F410M and F444W. At the time of writing, the JADES program is
still ongoing, and the parts of observations (PID: 1180 and 1210) included in
this study cover a total sky area of $\approx$ 60 arcmin2, including regions
around the Hubble Ultra-Deep Field where additional 5 medium-band NIRCam
imaging observations, i.e., F182M, F210M, F430M, F460M and F480M, have been
obtained by the JWST Extragalactic Medium-band Survey (JEMS, Williams et al.
2023, PID: 1963).

We reduce all these NIRCam imaging observations using a consistent method,
which has been presented in detail in Rieke & the JADES Collaboration (2023).
In short, we process the raw images using the JWST Calibration Pipeline with
some custom corrections. Specifically, during the Stage 1 of the JWST pipeline
we mask and correct for the “snowballs” effect caused by charge deposition
arising from cosmic ray hits. During the Stage 2 of the JWST pipeline we
remove the $1/f$ noise and subtract 2D background from the images, and also
correct the “wisp” features in the short-wavelength (SW, i.e. filters bluer
than F277W) images. Afterwards, we tie the astrometry of individual exposures
of a given visit to the World Coordinate System (WCS) of a reference catalog
constructed from HST/WFC3 F160W mosaics in GOODS-South with astrometry tied to
Gaia-EDR3 (Gaia Collaboration et al., 2021). Finally, we combine individual
visit-level mosaics together to create the final mosaic, where we choose a
pixel scale of 0.03 arcsec pixel-1 for all filters.

We also use the ancillary HST imaging observations taken from the latest
Hubble Legacy Fields Data Release in GOODS-South (Illingworth et al., 2016;
Whitaker et al., 2019), including the imaging data of the ACS/WFC F435W,
F606W, F775W, F850LP and F814W filters, and of the WFC3 F105W, F125W, F140W
and F160W filters. All these HST images are on the same WCS grid as our NIRCam
images.

##  3 Sample Selection

As the JADES program is still ongoing, the motivation of this work is not to
test, improve or develop selection methods to build a complete demographic
picture of quiescent galaxy populations across a wide dynamic range of stellar
mass. Instead, we focus this paper only on the population of massive quiescent
galaxies with stellar masses $M_{*}>10^{10}M_{\sun}$. Because existing HST
surveys are remarkably complete at this stellar-mass range ($>90\%$ complete,
Guo et al. 2013; Ji et al. 2018), our strategy is to make full use of the
power of those legacy surveys in GOODS-South to first select candidate
galaxies using methods that have been extensively tested in literature, and
then refine the selected galaxies using SED fitting with our new JWST data.

###  3.1 Selecting Candidate Quiescent Galaxies using Rest-UVJ Colors from
3D-HST

Our primary selection of quiescent galaxies is based on the redshift-dependent
criteria of rest-frame UVJ colors (Williams et al., 2009). We take advantage
of the measurements from the 3D-HST survey (Brammer et al., 2012; Skelton et
al., 2014). We first impose a cut on the signal-to-noise ratio (S/N $>10$) in
the HST/F160W filter, which ensures high-quality photometry and hence those
previous measurements of galaxy physical properties from 3D-HST. We then
select galaxies with $M_{*}>10^{9.7}M_{\sun}$111The SED fitting code
Prospector is elected to use in this work (Section 4.1). It has been noted
that the stellar mass measures from Prospector are $0.1-0.3$ dex higher (e.g.,
Leja et al., 2022; Ji & Giavalisco, 2022), therefore we use this
$10^{9.7}M_{\sun}$ cut to ensure our initial selection includes the vast
majority of galaxies with $M_{*}\geq 10^{10}M_{\sun}$ from Prospector. from
the 3D-HST catalog. Finally, we use the rest-frame $(U-V)$ and $(V-J)$ colors
from the 3D-HST catalog to make the initial selection of candidate quiescent
galaxies.

We select candidate quiescent galaxies whose rest-frame colors satisfy the UVJ
criteria from Schreiber et al. (2015):

$\begin{split}(U-V)>1.3,\\\ (V-J)<1.6,\\\ (U-V)>0.88(V-J)+0.49.\end{split}$
(1)

We add to our analysis additional seven222The most galaxies from Carnall et
al. (2020) have already been selected by the UVJ technique described above.
$z>2$ quiescent galaxies from Carnall et al. (2020) which appear in the JADES
footprint. Carnall et al. (2020) made the selection of $M_{*}>10^{10}M_{\sun}$
quiescent galaxies at $z>2$ based on their refined SED modelling that allows
sufficiently flexible solutions over a wide range of observed redshifts and
specific star formation rates (see their Section 2.3 and Table 2). Carnall et
al. (2020) showed that their selection is highly consistent with the UVJ
technique, with a 80% overlapping fraction, but is more complete in that their
SED modelling enabled the selection of relatively young quiescent galaxies
which have bluer $(U-V)$ colors that would otherwise have been missed using
the lower-redshift UVJ selection criteria.

The above selection gives us a total of 268 candidate massive quiescent
galaxies at $z>0.5$. We then remove galaxies hosting bright AGN from our
analysis, especially because the presence of AGN can affect the interpretation
of their morphological measures. We first remove X-ray AGN by cross matching
the sample with the X-ray catalog from the 7Ms Chandra Deep Field South (Luo
et al., 2017). We then remove IR AGN using IRAC colors (Donley et al., 2012;
Ji et al., 2022). We also remove radio AGNs identified by the latest and
deepest VLA radio observations in GOODS-South (Rujopakarn et al. 2023 in
preparation). Finally, we remove other types of AGN identified by recent SED
fitting from Lyu et al. (2022). Adding all these different AGN types together,
we identify a total of 71 AGN-hosting galaxies, corresponding to a fraction of
$26\pm 4\%$ (71/268) which is in broad agreement with that found in galaxies
of similar masses at Cosmic Noon (e.g. Xue et al., 2010; Ji et al., 2022). We
mention that we retain the $z_{\rm{spec}}=4.658$ quiescent galaxy GS-9209 that
has recently been spectroscopically confirmed with the JWST/NIRSpec
observation by Carnall et al. (2023b). This galaxy is not identified as an AGN
based on the selection methods considered above, but it does host a faint AGN,
as suggested by the broad H$\alpha$ emission in the NIRSpec spectrum. However,
because this AGN is so faint that its broad-band size measurement – except for
the morphology in F356W which covers H$\alpha$ (Section 6.2.2) – is not
significantly biased by AGN contribution, we decide to keep this galaxy in the
following analysis.

After removing the AGN-hosting galaxies, our initial selection based on 3D-HST
measurements yields a sample of 197 candidate quiescent galaxies at $0.5<z<5$.

###  3.2 Finalizing the Quiescent Sample using JADES Data

We use the new photometry from JADES to refine the initial sample selected
above. We perform SED fitting with a total of more than 20-filter photometry
from HST through JWST using the code Prospector (Johnson et al. 2021, see
Section 4.1 for details). We then use these newly obtained best-fit SEDs to
re-measure the rest-frame UVJ colors. At $z<2$, we retain the galaxies still
satisfying all three UVJ color cuts (Equation 1). At higher redshifts $z>2$,
quiescent galaxies are younger and have rest-frame UVJ colors more like post-
starburst galaxies (Whitaker et al., 2013; Belli et al., 2019; Carnall et al.,
2020). Following the finding of Carnall et al. (2020, see their Section 3.1),
we relax the horizontal cut on $(U-V)$ color and retain the galaxies
satisfying the remaining two UVJ color cuts. These give us a final sample of
161 robust UVJ-selected quiescent galaxies with $M_{*}>10^{10}M_{\sun}$ at
$0.5<z<5$, among which 130 and 31 galaxies are at $z<2$ and $z\geq 2$,
respectively.

Recent studies also used alternative selection methods of quiescent galaxies
based on a fixed or redshift-dependent threshold of specific star formation
rate (sSFR) (e.g., Speagle et al., 2014; Donnari et al., 2019; Leja et al.,
2022; Tacchella et al., 2022). However, accurately measuring sSFRs for
quiescent galaxies with photometry only is challenging (Conroy, 2013), which
is highly sensitive to SED model assumptions (e.g., Pacifici et al., 2023). As
Appendix A shows, we observe a significant systematic offset in the sSFR
measures by simply changing SFHs (Section 4.1), meaning that whatever sSFR
thresholds we use to select quiescent galaxies would be highly model-
dependent. By contrast, the rest-frame UVJ color measures are largely
insensitive to the SED model assumptions, making it a more empirical way to
select quiescent galaxies. Because of this, we choose to still use the UVJ
technique in this work. Nevertheless, the UVJ selection should be highly
consistent with the selection methods based on sSFR. To show this, in the left
panel of Figure 1 we plot the distribution of our final quiescent sample in
the rest-frame UVJ diagram. We color code each one of the galaxies according
to $(\frac{1}{\rm{sSFR}})/t_{\rm{H}}$, i.e., the ratio of the mass doubling
time at the current star formation rate333Here the current star formation rate
is defined as the star formation rate of the first lookback time bin, i.e.
0$-$30 Myr, from our fiducial SED fitting (Section 4.1). to the Hubble time at
the time of observation. A larger $(\frac{1}{\rm{sSFR}})/t_{\rm{H}}$
corresponds to a lower star-formation intensity. Regardless of SFH assumptions
(Appendix A), we find that the UVJ-selected quiescent galaxies have longer
mass doubling time than the Hubble time, confirming that the UVJ-selected
quiescent galaxies also have low sSFRs.

Figure 1: The final sample of 161 massive quiescent galaxies presented in this
work (Section 3). Quiescent galaxies at $z>2$ are shown as stars, while those
at $0.5<z<2$ are shown as circles. The orange square marks the quiescent
galaxy GS-9209 at $z_{\rm{spec}}=4.658$ (Carnall et al., 2023b). Each one of
the galaxies is color-coded according to its star formation intensity,
quantified using $(\frac{1}{\rm{sSFR}})/t_{\rm{H}}$ (Section 3.2). Left: Rest-
frame UVJ color diagram. Background grey contours show the distribution of all
galaxies with stellar mass $\log M_{*}>10$ and redshift $0.5<z<5$ from 3D-HST.
Black dashed lines show the criteria from Schreiber et al. (2015), which we
use for the selection of quiescent galaxies at $z<2$. Following recent
studies, we relax the criterion of the horizontal cut on $(U-V)$ to select
quiescent galaxies at $z\geq 2$ (Section 3). Right: The distribution of the
quiescent galaxy sample in the stellar mass-redshift space.

##  4 Measurements

###  4.1 SED Fitting with Prospector

With the JADES photometry that samples the wavelength range from rest-frame UV
to NIR with more than 20 filters, we perform SED fitting using the code
Prospector (Johnson et al., 2021), with the emphasis on robustly measuring
stellar mass, photometric redshift and sSFR. We assume a similar SED model
used in Ji et al. (2023). In brief, we adopt the FSPS stellar synthesis code
(Conroy et al., 2009; Conroy & Gunn, 2010) with the stellar isochrone
libraries MIST (Choi et al., 2016; Dotter, 2016) and the stellar spectral
libraries MILES (Falcón-Barroso et al., 2011). We adopt the Madau (1995) IGM
absorption model and include the nebular emission model of Byler et al.
(2017). We model the dust attenuation following Tacchella et al. (2022) where
the dust attenuation of nebular emission and young stellar populations, and of
old stellar populations are treated differently (Charlot & Fall, 2000). We fix
the redshift when a spectroscopic one is available444We use the spectroscopic
redshift catalog in GOODS-South compiled by the ASTRODEEP project (Merlin et
al., 2021)., otherwise we fit it as a free parameter with a flat prior of
$z\in(0.1,10)$.

As the fiducial model of this study, we assume nonparametric, piece-wise SFH
composed of 9 lookback time bins, where SFR is constant within each bin. We
model nonparametric SFH using the Dirichlet prior (Leja et al., 2017) which
can recover high-fidelity SFHs of different galaxy types while effectively
mitigate the overfitting problem (Leja et al., 2019). In addition, this prior
has been demonstrated, with synthetic observations from cosmological
simulations, to give the unbiased inference of stellar ages for massive
quiescent galaxies at $z\sim 2$ (see Appendix A in Ji & Giavalisco, 2022).

In Figure 2 we compare the redshifts and stellar masses from our Prospector
fitting with those from 3D-HST (Brammer et al., 2012; Skelton et al., 2014).
The agreement between the redshift measures is excellent. Our Prospector
measurement of stellar masses is also tightly correlated with that from
3D-HST, although our stellar masses are systematically higher by $\approx
0.16$ dex, a quantitatively similar systematics has also been found by other
studies (e.g., Leja et al., 2019; Lower et al., 2020; Leja et al., 2022; Ji &
Giavalisco, 2022, 2023). In Appendix A, we further test our measurements with
another two different SFHs, namely nonparametric SFHs with the continuity
prior (Leja et al., 2019) and delayed-tau parametric SFHs. We show that our
measurements of stellar mass, rest-frame UVJ colors (hence the selection of
quiescent galaxies) are robust.

Figure 2: Comparisons of the redshift (top) and stellar-mass (bottom)
measures. Galaxies with spectroscopically confirmed redshifts are color coded
in pink. GS-9209 is marked with the orange square. The grey, dashed line marks
the one-to-one relationship. The uncertainty of photometric redshifts from our
Prospector fitting is smaller than the point, hence hard to see in the plot.
Our stellar-mass measures with Prospector are systematically higher by 0.16
dex relative to the measures from 3D-HST.

###  4.2 Measuring Rest-frame Sizes

In this paper we refer to the circularized effective (half-light) radius
($R_{e}=R_{e,maj}\times\sqrt{b/a}$) as the size of galaxies555We find no
substantial changes in our conclusions if we use the effective semi-major axis
instead.. We measure the size at three characteristic rest-frame wavelengths –
0.3, 0.5 and 1$\micron$ – using different wide-band images. As Figure 3
illustrates, we can perform the most measurements using NIRCam images from
JADES, except the rest-0.3$\micron$ sizes for galaxies at $z<1.3$ which we
instead use the HST/ACS F606W images. We have tested our results by linearly
interpolating666For the two $z>4$ galaxies in our sample, we actually need do
the extrapolation to get their rest-frame $1\micron$ sizes. the measurements
of two adjacent filters to get the size at the exact rest-frame wavelength. We
found no substantial changes in any of our conclusions.

Figure 3: Wide-band filters used for the rest-frame 0.3$\micron$ (blue),
0.5$\micron$ (green) and 1$\micron$ (red) size measures. Expect the HST/ACS
F606W, all other filters are from JWST/NIRCam.

We measure the size by performing PSF-convolved morphological fitting using
the software Galfit (version 3.0.5, Peng et al., 2010a), where each one of the
galaxies is modelled as a single 2-dimensional Sérsic profile. Before running
Galfit, we first center on each one of the galaxies to make a
$7^{\prime\prime}\times 7^{\prime\prime}$ cutout. Similar to previous studies
(e.g., van der Wel et al., 2012), we only allow the Sérsic index $n$ varying
between 0.2 and 8. To mitigate possible impacts from adjacent galaxies on the
morphological fitting, we use the JADES segmentation map (Rieke & the JADES
Collaboration, 2023) to mask out all neighboring galaxies in the cutouts. By
default, we fit the local sky background of each cutout as a free parameter.
The PSFs used for the fitting are the model PSFs (mPSFs) of JADES (Ji et al.,
2023). The mPSFs are generated using the software Webbpsf (v1.1.1, Perrin et
al., 2014a), where we have carefully taken into account the effect of data
reduction. As demonstrated by Ji et al. (2023, see their Appendix A), the
mPSFs are highly consistent with the commonly-used effective PSFs (ePSFs,
Anderson & King, 2000) to the level of $\approx 1\%$ in terms of enclosed
energy radial profiles. The biggest advantage of using mPSFs, especially for
the extragalactic deep fields like JADES which contain very limited number of
bright, unsaturated isolated stars in their footprints, is that it allows us
to have robust, stable PSF models to larger ($>1^{\prime\prime}$) angular
scales.

We estimate uncertainties of the fitted parameters using a Monte Carlo method.
Specifically, we use the error flux extension of the NIRCam data to Monte
Carlo resample the image pixel values 300 times. We then run Galfit on the
resampled images and use the range between 16th and 84th percentiles as
$1\sigma$ uncertainties. In this paper, because we only focus on massive
galaxies which are detected with S/N $>$ 100 in the JWST images777For the
galaxies presented in this paper, the median S/N is 220 in F150W, and is 450
in F444W, where S/N is measured within the Kron aperture., the derived
uncertainty of $R_{e}$ is generally very small. For reference, we find a
median uncertainty of NIRCam/F150W size to be $\delta R_{e}/R_{e}\approx 2\%$,
which is in quantitative agreement with the uncertainties estimated for
galaxies with similar S/Ns using HST/F160W images (see Figure 6 in van der Wel
et al. 2012 and Figure 6 in Ji et al. 2022).

We now test systematics of our Galfit fitting. First, we consider the
systematics introduced by the local sky background estimation. Instead of
fitting the sky as a free parameter, if we fix it to the $3\sigma$-clipped
median pixel value of the cutout after masking out all detected sources using
the segmentation map, the median relative difference in the size measures will
be $\approx 1\%$. Second, instead of using mPSFs, if we use ePSFs in the
Galfit fitting, the median relative difference in the size measures will be
$\approx 3\%$. These systematic uncertainties are small and have no impacts on
the results of this work. Finally, we compare our NIRCam/F150W sizes to the
HST/WFC3 F160W sizes measured by van der Wel et al. (2012). As Figure 4 shows,
we find excellent agreement between the two measures, with a median difference
of $0.03\pm 0.07$ dex. Moreover, among the 161 quiescent galaxies, we find
that 24 F150W Galfit fits hit the Sérsic index bounds ($n=0.2$ or 8), which is
also consistent with van der Wel et al. (2012) who found 21 F160W fits having
the same issue.

We remove galaxies with $n=0.2$ or 8 from our sample to ensure reliable size
measures for the size evolution analysis. Specifically, we remove 37, 25 and
25 galaxies from the analysis of the rest-frame 0.3$\micron$, 0.5$\micron$ and
1$\micron$ size evolution, respectively. We note that all of the removed fits
have $n=8$, indicating that those galaxies have very compact cores. Removing
them from our analysis thus can drag the median size up. To check the
magnitude of this impact, for the galaxies with unreliable single Sérsic fits,
we use the alternative approach with the Richardson–Lucy deconvolution
algorithm (see Section 4.2.1 below for details) to estimate their sizes, and
then add these alternative sizes back to the analysis. Despite that the best-
fit parameters of the size evolution are consistent with those reported in
Table 1 within uncertainty, the intercept $C$ (Equation 3) systematically
reduces by 0.08, 0.06 and 0.05 dex for rest-frame 0.3$\micron$, 0.5$\micron$
and 1$\micron$ sizes, respectively. Other than this, we find no substantial
changes in any of our results.

Figure 4: Size from our measurements using the NIRCam/F150W images ($y$-axis,
Section 4.2) as a function of size from van der Wel et al. (2012) who
performed the measurements using the HST/WFC3 F160W images ($x$-axis). The
agreement between the two measurements is great, with a median difference
($y-x$) of 0.03 dex and a standard deviation of 0.07 dex.

####  4.2.1 Testing the size measures for the $z>3$ galaxies

As we will show in Section 5.2, the $z>3$ quiescent galaxies in our sample
have typical sizes of $R_{e}^{1\micron}\sim 0.3$ kpc, corresponding to 0”.04
at $z=3$ which is only slightly larger than the 0”.03 pixel scale of JADES
mosaics. Moreover, all the $z>3$ quiescent galaxies have large Sérsic indices
$n>2$. Long tails of large Sérsic index profiles can make the analysis
sensitive to profile mismatch and sky background estimates (Peng et al.,
2010a). We thus conduct several further tests to check the robustness of the
rest-frame 1$\micron$ size estimates for these extremely compact galaxies at
$z>3$.

First, the robustness of Sérsic fitting, especially for galaxies with very
compact morphologies, relies on accurate pixelization and integration of
Sérsic profiles. As shown by Robotham et al. (2017, see their Figure 7), the
flux weighted error model image generated by Galfit increases with Sérsic $n$.
We thus test our Sérsic fitting for the $z>3$ quiescent galaxies with
Lenstronomy (Suess et al., in preparation), a newly developed morphological
fitting code that adopts new algorithms to generate Sérsic models and to
conduct fitting. We plot the results as black circles in Figure 5. Excellent
agreement is seen between the sizes from Galfit and from Lenstronomy.

Second, the degeneracy between $R_{e}$ and Sérsic $n$ (e.g. Ji et al., 2020,
2022) makes the size measures uncertain. To check how sensitive the size
estimates are due to uncertain $n$, instead of fitting $n$ as a free
parameter, we run Galfit by fixing $n=1$ (exponential disk) or $n=4$ (de
Vaucouleurs profile) for the $z>3$ quiescent galaxies. In Figure 5 the results
are shown as blue and red squares, respectively. The sizes from the $n=4$
Sérsic fitting are in very good agreement with our fiducial measurements,
which is expected since all the galaxies have $n>2$ from the default Sérsic
fitting. The agreement between the fiducial sizes and those from the $n=1$
fitting is also good, with the differences well within a factor of 2. We note
that, relative to our default Sérsic fitting, the reduced $\chi^{2}$ increases
by a factor of 1.8 and 1.1 for the fitting with fixed $n=1$ and $n=4$,
respectively, showing that the quality of fitting is better with $n$ as a free
parameter, as expected.

Third, it is known that the light profile of quiescent galaxies can be more
complex than a single Sérsic profile (e.g Huang et al., 2013a, b). Davari et
al. (2014) simulated quiescent galaxies at different redshifts by mimicing
their light profiles using the detailed substructures observed in nearby
elliptical galaxies. They found that the systematic uncertainty of the size
measures at $z\sim 2$ from single Sérsic fitting is at the level of $\sim
10-20$%. If high-redshift quiescent galaxies indeed have similarly multi-
component light profiles like lower-redshift counterparts, the tests done by
Davari et al. (2014) will then suggest that the mismatch of assumed light
profiles only introduces subtle bias in the inferred sizes. However, because
the rest-frame NIR light distribution of quiescent galaxies at $z>3$ remains
largely unknown, we further test our size measures using the following
nonparametric approach with the Richardson–Lucy deconvolution algorithm
(Richardson, 1972).

Image deconvolution is mathematically ill-posed. Yet, with precise PSF models,
images can be restored to some accuracy through an iterative Richardson–Lucy
algorithm. As a sanity check, we run Richardson–Lucy deconvolution for the
mPSF image (Section 4.2) of F444W, i.e. the filter used for the
rest-1$\micron$ size measures at $z>3$ (Figure 3). After $\sim 50$ iterations,
the FWHM of the deconvolved mPSF image asymptotes to 0”.015, i.e. half the
pixel size which is the theoretical limit we can achieve for JADES images with
the Richardson–Lucy deconvolution, showing the effectiveness of this
technique. Thus our strategy is to first deconvolve the F444W images of the
$z>3$ quiescent galaxies, and then directly estimate the effective radii
($R_{e}^{\rm{RL}\dagger}$) of the galaxies by measuring their radial light
profiles in the deconvolved images. In this way the PSF effect can be largely
mitigated, and the difference between $R_{e}^{\rm{RL}\dagger}$ of the galaxies
and that of the PSF, i.e. 0”.015, can be considered as a result of the
intrinsic broadening due to the sizes of the galaxies. The sizes can then be
inferred following

$R_{e}^{\rm{RL}}[\rm{arcsec}]=\sqrt{(R_{e}^{\rm{RL}\dagger})^{2}-(0.015)^{2}}$
(2)

The biggest advantage of the method above based on the Richardson–Lucy
deconvolution is that it does not require an assumed form for the parametric
light profile, and any complex structures in the galaxies are naturally taken
into account during the $R_{e}^{\rm{RL}}$ estimates. Moreover, this method
does not bias towards brighter parts of the galaxies (which differs from
Galfit which conducts $\chi^{2}$ minimization fitting and hence essentially
provides light-weighted measures). The disadvantage of this method, though, is
that no formal stop criterion exists for the iterative process (Prato et al.,
2012). Previous HST studies by van der Wel et al. (2011) fixed the number of
iterations to be 16. In our cases, we find that, once the iteration number is
greater than 10, the $R_{e}^{\rm{RL}}$ measures become very stable for the
$z>3$ quiescent galaxies.

In Figure 5, $R_{e}^{\rm{RL}}$ are shown as orange triangles (based on 50
times Richardson–Lucy iterations). Even with this very different approach, we
still find very good agreement with the sizes derived from the default single
Sérsic fitting with Galfit. We also note that, despite the $z>3$ quiescent
galaxies are very compact, our analysis based on the Richardson–Lucy
deconvolution suggests that their light profiles are inconsistent with point
sources, as all of them have $R_{e}^{\rm{RL}\dagger}>0".015$, i.e. the
galaxies are (marginally) resolved in F444W.

Figure 5: Rest-frame 1$\micron$ sizes of the $z>3$ quiescent galaxies from
different methods. The $x$-axis shows the fiducial measurements in this work
(Section 4.2). The $y$-axis shows the measurements from alternative methods
detailed in Section 4.2.1, including using (1) a different morphological
fitting tool Lenstronomy (black circles), (2) Galfit but fixing Sérsic index
$n=1$ (blue squares) or $n=4$ (red squares) during the fit and (3) a
nonparametric method with the Richardson–Lucy deconvolution (orange
triangles). The main purpose here is to compare the sizes from different
methodologies, we thus do not estimate uncertainties for the $y$-axis. The
dashed line marks the one-to-one relation, and the dotted lines mark the 2
times above/below the one-to-one relation. For the vast majority of the $z>3$
quiescent galaxies, the relative difference in sizes from different methods is
$<50\%$.

Finally, we stress that the purpose of running different tests here is not to
evaluate which methodology is the best in characterizing the morphology of
high-redshift quiescent galaxies. Instead, our goal is to check if the size
estimates of the $z>3$ quiescent galaxies suffer from any significant
systematic uncertainties from the imperfect fitting algorithm, the degeneracy
between $R_{e}$ and $n$, and the mismatch of assumed light profiles.
Regardless of methods, we find strong correlations between the sizes from our
fiducial measurements and the sizes from other alternative methods described
above. The difference in sizes from different methods are well within a factor
of 2. In fact, the relative difference in size is $<50\%$ for the vast
majority of the $z>3$ galaxies (Figure 5). We thus conclude that the size
estimates for the quiescent galaxies at $z>3$ in our sample are not sensitive
to significant systematic uncertainties, hence the strong size evolution of
quiescent galaxies towards $z>3$ reported below (Section 5.2) is robust.

##  5 Results

Figure 6 presents the RGB images of individual quiescent galaxies in this
work, where the RGB filters are chosen to probe rest-frame 0.3, 0.5, and
1.0$\micron$ (Figure 3) and each galaxy is plotted at the corresponding
location in the diagram of rest-1$\micron$ size vs redshift. We observe strong
size evolution of quiescent galaxies from $z\sim 5$ to $z\sim 0.5$. In what
follows, we will quantify and compare the size evolution at different rest-
frame wavelengths from UV through NIR. We will also explore the possible
dependence of size evolution on the physical properties of quiescent galaxies.
Finally, we will present the stellar mass-size relationship of massive
quiescent galaxies.

Figure 6: RGB cutouts of massive quiescent galaxies presented in this work.
Each one of the galaxies is plotted at its location in the diagram of rest-
frame 1$\micron$ size evolution. The filters used to produce the RGB cutouts
are illustrated in Figure 3. The size of the cutouts is 20 kpc $\times$ 20
kpc.

###  5.1 Multi-wavelength Sizes: Piecing Together the Physical Drivers of Size
Evolution

Because of the spatial variations of dust attenuation and stellar-population
properties, galaxies show radial gradients of stellar mass-to-light ratio,
which makes the interpretation of the apparent evolution of projected light
profiles nontrivial (e.g. Suess et al., 2019a, b; Mosleh et al., 2020; Miller
et al., 2022; van der Wel et al., 2023). In this work, we perform the
morphological analysis at three different wavelengths – that are sensitive to
different physical processes – to constrain the physical drivers of the size
evolution of quiescent galaxies. The two shorter wavelengths are chosen to be
rest-frame 0.3$\micron$ and 0.5$\micron$, which cover key stellar population
features around the rest-frame 4000Å region. The third characteristic
wavelength is chosen to be rest-frame 1$\micron$ which is a sensitive probe of
stellar-mass distributions, and is the longest wavelength where we can still
measure galaxy morphology with the JADES NIRCam imaging over the redshift
range considered here (Figure 3).

Compared to rest-UV and optical, rest-NIR light profiles are a much better
proxy for the stellar-mass distribution of galaxies, because the stellar mass-
to-light ratio at NIR only weakly, if at all, depends on recent star formation
(e.g., Bell & de Jong, 2001; McGaugh & Schombert, 2014). To show this, we
cross match our quiescent sample with the catalog from Suess et al. (2020) who
measured galaxy half-mass sizes ($R_{\rm{halfmass}}$) using HST data based on
spatially resolved SED fitting. A total of 116 galaxies in our sample have the
$R_{\rm{halfmass}}$ measures. In the top panel of Figure 7, we plot
$R_{\rm{halfmass}}$ as a function of $R_{e}$ measured at different
wavelengths. Strong correlations are observed between $R_{\rm{halfmass}}$ and
$R_{e}$, regardless of which wavelength $R_{e}$ is measured. Systematic
offsets are also observed between $R_{\rm{halfmass}}$ and $R_{e}$, likely due
to the different methodologies used for the size measures.

Figure 7: Half-light sizes $R_{e}$ vs half-mass sizes $R_{\rm{halfmass}}$ for
the subsample of quiescent galaxies whose $R_{\rm{halfmass}}$ measures are
available from Suess et al. (2019b). Rest-frame 0.3, 0.5 and 1$\micron$
$R_{e}$ are shown in blue, green and red, respectively. The bottom panel shows
histograms of the difference between the half-light and half-mass sizes
($\Delta$), and the inset shows the standard deviation of $\Delta$ (Section
5.2). The agreement between $R_{\rm{halfmass}}$ and rest-1$\micron$ $R_{e}$ is
the best.

We compare the agreement between $R_{\rm{halfmass}}$ and $R_{e}$ measured at
different wavelengths. In the bottom panel of Figure 7, we show the histograms
of $\Delta$, the relative difference between $R_{\rm{halfmass}}$ and $R_{e}$.
We measure the standard deviation of $\Delta$, i.e. $\sigma(\Delta)$. We
estimate the uncertainty of $\sigma(\Delta)$ through a bootstrap Monte Carlo
method. Specifically, we bootstrap the sample, use normal distributions to
Monte Carlo resample the values of $R_{\rm{halfmass}}$ and $R_{e}$ with the
corresponding measurement uncertainties, and finally we measure
$\sigma(\Delta)$. We repeat the procedure 1000 times and use the range between
16th and 84th percentiles as the uncertainty of $\sigma(\Delta)$. As the inset
of the bottom panel of Figure 7 shows, the $\Delta$ between
$R_{\rm{halfmass}}$ and $R_{e}^{1\micron}$ has the smallest scatter, showing
rest-NIR light profiles indeed a robust, much better proxy for stellar-mass
distributions than the rest-UV and optical ones.

In what follows, we thus will use $R_{e}^{1\micron}$ as the size proxy for
stellar-mass distributions of the galaxies. Because the 4000Å break contains
key information related to dust attenuation, stellar age and metallicity,
combining with the size evolution at rest-frame 0.3$\micron$ and 0.5$\micron$,
we will be able to empirically constrain the physical mechanisms driving the
size evolution of quiescent galaxies.

###  5.2 Size Evolution of Quiescent Galaxies from UV through NIR

We characterize the size evolution by fitting the relationship between the
$R_{e}$ and $z$ of individual quiescent galaxies with a functional form of

$\log R_{e}=\beta\log(1+z)+C.$ (3)

To estimate the uncertainties of $\beta$ and $C$, we bootstrap the entire
sample, use normal distributions to Monte Carlo resample the size and redshift
measures with the corresponding uncertainties, and finally fit the size
evolution. We repeat this bootstrap Monte Carlo procedure 1000 times and use
the range between 16th and 84th percentiles as $1\sigma$ uncertainties. The
best-fit relationships are present in Table 1 and Figure 8.

Figure 8: Rest-frame size evolution of massive quiescent galaxies. The
corresponding best-fit parameters can be found in Table 1. Top: The left panel
shows the best-fit size evolution at rest-frame 0.3$\micron$ (blue),
0.5$\micron$ (green) and 1$\micron$ (red), respectively. The shaded regions
mark the $1\sigma$ uncertainties of the best-fit relations. The filled squares
with error bars are median $R_{e}$ and their uncertainties in individual
redshift bins (Section 5.2). The right panel shows the 1-, 2- and 3-$\sigma$
contours of the $R_{e}$-z relations based on our bootstrap Monte Carlo method
(Section 5.2). Bottom: Differing from the top-left panel, here we plot
individual quiescent galaxies to the size-evolution diagram. Like previous
figures, GS-9209 is marked with the orange square.

####  5.2.1 The pace of the size evolution

Regardless of rest-frame wavelengths, the size of quiescent galaxies evolves
with redshift at a similar pace, i.e., a similar $\beta\approx-1.25\pm 0.20$
(Table 1). Because the fitting of size evolution was done using the
measurements of individual galaxies and most of the quiescent galaxies in our
sample are at $z<2$, the low-redshift quiescent galaxies have more weight in
the fitting procedure. Notwithstanding the very limited sample size at $z>3$,
by looking at the bottom panels of Figure 8 we note that the vast majority of
$z>3$ quiescent galaxies are below (i.e. have smaller sizes than) the best-fit
$R_{e}$-z relationship. To better show this, we estimate the median sizes of
quiescent galaxies at different redshifts. We divide the entire sample into
five redshift bins, i.e., four bins between $z=0.5$ and $z=2.5$ with an
interval of $\Delta z=0.5$ and the last bin of $z>2.5$. Similarly, we use the
bootstrap Monte Carlo method to estimate the uncertainties of the median
sizes. As the top-left panel of Figure 8 shows, the median size of the highest
redshift bins deviates from the best-fit size evolution – quiescent galaxies
at $z>3$ are smaller than the predictions by extrapolating the size evolution
of lower-z quiescent galaxies. We will discuss the implications regarding this
finding later in Section 6.2.

Table 1: Median sizes of massive quiescent galaxies and their best-fit evolution with redshift. | $z=0.5\sim 1.0$ | $z=1.0\sim 1.5$ | $z=1.5\sim 2.0$ | $z=2.0\sim 2.5$ | $z=2.5\sim 5.0$ |  $\beta$(b) |  $C$(b)
---|---|---|---|---|---|---|---
$R_{e}(0.3\micron)^{(a)}$ | 2.11${}_{-0.33}^{+0.52}$ | 1.39${}_{-0.22}^{+0.26}$ | 1.47${}_{-0.31}^{+0.21}$ | 0.94${}_{-0.08}^{+0.43}$ | 0.56${}_{-0.07}^{+0.19}$ | -1.23${}_{-0.24}^{+0.21}$ | 0.62${}_{-0.11}^{+0.10}$
$R_{e}(0.5\micron)^{(a)}$ | 1.90${}_{-0.26}^{+0.17}$ | 1.30${}_{-0.13}^{+0.16}$ | 1.07${}_{-0.14}^{+0.17}$ | 0.72${}_{-0.06}^{+0.42}$ | 0.60${}_{-0.21}^{+0.11}$ | -1.24${}_{-0.17}^{+0.21}$ | 0.57${}_{-0.08}^{+0.07}$
$R_{e}(1.0\micron)^{(a)}$ | 1.55${}_{-0.20}^{+0.17}$ | 1.10${}_{-0.13}^{+0.13}$ | 0.96${}_{-0.09}^{+0.19}$ | 0.6${}_{-0.05}^{+0.07}$ | 0.38${}_{-0.04}^{+0.17}$ | -1.33${}_{-0.18}^{+0.19}$ | 0.53${}_{-0.07}^{+0.06}$
$\dagger$ Values shown in the table are the medians and their 1-$\sigma$
uncertainties from our bootstrap Monte Carlo method (Section 5.2).
(a) In the unit of kpc. (b) The best-fit relation of $\log
R_{e}=\beta\log(1+z)+C$.

We continue to divide the sample into two mass bins using the median stellar
mass $\langle\log(M_{*}/M_{\sun})\rangle=10.6$. We present the size evolution
of the two mass bins in Figure 9 and Table 2. We observe a dependence of the
size evolution on stellar mass, where the evolution is faster – more negative
$\beta$ (Equation 3) – for more massive quiescent galaxies. This finding does
not depend on the rest-frame wavelength at which the sizes are measured (Table
2).

We further check the robustness of this stellar-mass dependence. Considering
the relatively small sample size of this study (161 galaxies), instead of
dividing the sample into more stellar mass bins, we first sort stellar masses
of individual quiescent galaxies into an increasing order. Then, starting from
the first 30 of the sorted sample, we keep adding more massive quiescent
galaxies into the fit of size evolution. We finally study the change of
$\beta$ as a function of the maximum stellar mass of quiescent galaxies
included in the fit. As Figure 10 shows, $\beta$ becomes increasingly negative
as more massive quiescent galaxies included, confirming that the size
evolution is faster for more massive quiescent galaxies. In Section 6.1, we
will discuss in detail the implications of this finding for the physical
mechanisms governing the size evolution of massive quiescent galaxies.

Figure 9: Stellar-mass dependence of the size evolution at rest-1$\micron$. The left panel shows the size evolution of quiescent galaxies with $\log(M_{*}/M_{\sun})=10-10.6$, while the right shows the evolution for those with $\log(M_{*}/M_{\sun})>10.6$. Little dots are individual quiescent galaxies, and large squares are the median size of each redshift bin (the same as Figure 8). The dashed red line shows the best-fit evolution (Equation 3), and the shaded region marks the 1$\sigma$ uncertainty estimated using the bootstrap Monte Carlo method. The size evolution is faster for more massive quiescent galaxies. The similar is also found for the size evolution at rest-0.3$\micron$ and 0.5$\micron$ (Table 2). Like previous figures, GS-9209 is marked with the orange square.  Figure 10: The change of the best-fit $\beta$ (Equation 3) as more massive quiescent galaxies are included in the fit (black solid line). The x-axis shows the maximum stellar mass of the quiescent galaxies included in the fit (see Section 5.2.1 for details). The gray shaded regions mark the 1$\sigma$ uncertainty estimated using the bootstrap Monte Carlo method. As more massive quiescent galaxies add to the fit, the slope $\beta$ becomes more negative, i.e. the size evolution is faster for more massive quiescent galaxies. Table 2: Stellar-mass dependence of the size evolution of massive quiescent galaxies (Equation 3). $\lambda_{\rm{rest}}$ | $M_{*}=10^{10}\sim 10^{10.6}M_{\sun}$ | $M_{*}>10^{10.6}M_{\sun}$
---|---|---
$\beta$ | $C$ | $\beta$ | $C$
$0.3\micron$ | -1.06${}_{-0.33}^{+0.31}$ | 0.50${}_{-0.12}^{+0.12}$ | -1.68${}_{-0.28}^{+0.27}$ | 0.87${}_{-0.14}^{+0.14}$
$0.5\micron$ | -1.04${}_{-0.26}^{+0.25}$ | 0.43${}_{-0.09}^{+0.08}$ | -1.60${}_{-0.22}^{+0.29}$ | 0.79${}_{-0.11}^{+0.09}$
$1\micron$ | -1.15${}_{-0.22}^{+0.25}$ | 0.38${}_{-0.07}^{+0.08}$ | -1.73${}_{-0.20}^{+0.24}$ | 0.76${}_{-0.09}^{+0.09}$

####  5.2.2 The difference in size at different wavelengths

The intercept $C$ (Equation 3) of the size evolution depends on rest-frame
wavelengths (Figure 8). Massive quiescent galaxies are larger at shorter
wavelengths: relative to their rest-1$\micron$ sizes, on average, quiescent
galaxies are larger by 45% at rest-frame 0.3$\micron$ and 15% at rest-frame
0.5$\micron$. Our results are in quantitative agreement with other recent JWST
studies (Suess et al., 2022; van der Wel et al., 2023), showing that the
stellar-mass distribution of quiescent galaxies is significantly more compact
than their rest-frame UV and optical light distributions, caused by the radial
changes in stellar-population properties such as stellar age, metallicity and
dust attenuation (e.g., Suess et al., 2022; Miller et al., 2023; van der Wel
et al., 2023; Suess et al., 2023).

We attempt to constrain the key physical quantity driving the differences in
morphology at different rest-wavelengths. In Figure 11, we plot the rest-NIR
to optical size ratio, i.e. R${}_{e}^{1\micron}$/R${}_{e}^{0.5\micron}$, as a
function of the mass-weighted stellar age. Because the strength of color
gradients depends on stellar mass (Tortora et al., 2010; Suess et al., 2019b;
van der Wel et al., 2023), we separate our sample into two stellar-mass bins
using $\langle\log(M_{*}/M_{\sun})\rangle=10.6$, the median value of the
entire sample. Regardless of stellar mass, younger quiescent galaxies with
ages $\leq 1.5$ Gyr have relatively flat color gradients, i.e.
R${}_{e}^{1\micron}$/R${}_{e}^{0.5\micron}$ closer to 1, in broad agreement
with studies of post-starburst galaxies (Maltby et al., 2018; Suess et al.,
2020). As galaxies become older ($>1.5$ Gyr), we observe an increase in the
number of quiescent galaxies with large difference in rest-optical and NIR
morphologies (R${}_{e}^{1\micron}$/R${}_{e}^{0.5\micron}<0.8$). Color
gradients of older quiescent galaxies are more negative, i.e. their centers
are redder. This is consistent with the general picture of inside-out growth
of massive galaxies, followed by inside-out quenching, where galaxies built up
their centers earlier than outskirts (e.g. Tacchella et al., 2015; Nelson et
al., 2016; Ji & Giavalisco, 2023).

We see tentative evidence that old ($>$1.5 Gyr), higher-mass quiescent
galaxies have flatter color gradients than lower-mass ones. The median size
ratio $\langle$R${}_{e}^{1\micron}$/R${}_{e}^{0.5\micron}$$\rangle$ is
0.87$\pm$0.02 for $\log(M_{*}/M_{\sun})>10.6$ quiescent galaxies, whereas it
is 0.83$\pm$0.02 for $10<\log(M_{*}/M_{\sun})<10.6$ ones (the quoted
uncertainty is estimated following the bootstrap Monte Carlo method). This is
in line with the studies of nearby massive ($\log(M_{*}/M_{\sun})>10$)
galaxies from SDSS: Tortora et al. (2010) found that the most massive
($\log(M_{*}/M_{\sun})>11$) early type galaxies have almost flat color
gradients which become increasingly more negative as early type galaxies
become less massive. This finding needs to be tested with future larger
samples of massive quiescent galaxies at high redshifts.

We also check the relationship of R${}_{e}^{1\micron}$/R${}_{e}^{0.5\micron}$
with dust attenuation $\rm{E(B-V)}$ and stellar metallicity $Z_{*}$. We do not
see clear trends, although this may be a result of the degeneracy between
$\rm{E(B-V)}$ and $Z_{*}$ that makes the inference of them individually highly
uncertain with photometric data only. We defer detailed analysis of the
spatially resolved stellar populations of the quiescent galaxies to a future
study aimed at dissecting the physical mechanism(s) behind this finding.

Finally, in addition to R${}_{e}^{1\micron}$/R${}_{e}^{0.5\micron}$, we have
also studied the trends discussed above using
R${}_{e}^{1\micron}$/R${}_{e}^{0.3\micron}$. We do not find any substantial
changes in our conclusions, although the measurement of
R${}_{e}^{1\micron}$/R${}_{e}^{0.3\micron}$ is noisier than
R${}_{e}^{1\micron}$/R${}_{e}^{0.5\micron}$ as quiescent galaxies are fainter
in rest-UV than in rest-optical.

Figure 11: The ratio of rest-1$\micron$ size to rest-0.5$\micron$ size as a
function of mass-weighted stellar age. In the top panel galaxies are color
coded according to the color excess $\rm{E(B-V)}$, while in the bottom panel
galaxies are color coded according to the stellar metallicity $Z_{*}$. The
horizontal dashed line marks R${}_{e}^{1\micron}$/R${}_{e}^{0.5\micron}=1$,
i.e. no evidence of color gradients. The black solid line marks
R${}_{e}^{1\micron}$/R${}_{e}^{0.5\micron}=0.8$, about the median value of all
quiescent galaxies in this study. Younger quiescent galaxies in general have
rest-optical to rest-NIR size ratios closer to unity, i.e. an indication of
flatter color gradients. As they become older, we observe an emergence of
quiescent galaxies with large color gradients, i.e.
R${}_{e}^{1\micron}$/R${}_{e}^{0.5\micron}<0.8$.

####  5.2.3 Comparing the size evolution of this work with other studies

We first compare our results with HST studies. van der Wel et al. (2014)
combined all five CANDELS fields to measure the rest-0.5$\micron$ size
evolution over $0<z<2.5$. Our new measurements obtain $\beta\approx-1.25\pm
0.20$ (Table 1) which is broadly consistent with van der Wel et al. (2014) who
found $-1.3\lesssim\beta\lesssim-1$ for $\log(M_{*}/M_{\sun})>10$ quiescent
galaxies (see their Table 2). This shows that the size evolution of quiescent
galaxies established based on rest-optical HST imaging still holds with the
deeper, higher resolution JWST/NIRCam observations.

More specifically, regarding the size evolution of quiescent galaxies with
$\log(M_{*}/M_{\sun})=10-10.6$, we obtain $\beta=-1.04^{+0.25}_{-0.26}$ (Table
2), which is in quantitative agreement with van der Wel et al. (2014) who
found $\beta=-1.01\pm 0.06$ for quiescent galaxies of similar stellar masses.
Regarding the higher-mass ones, i.e. $\log(M_{*}/M_{\sun})>10.6$, we obtain
$\beta=-1.60^{+0.29}_{-0.22}$, compared to $\beta=-1.32\pm 0.21$ from van der
Wel et al. (2014). Statistically, this difference in $\beta$ is marginal,
$\lesssim 1\sigma$. Nonetheless, we still investigate what causes the
marginally steeper size evolution we found for the $\log(M_{*}/M_{\sun})>10.6$
quiescent galaxies. First, this is not due to the systematics in the size
measures between JWST and HST, as we found consistent sizes between the
measurements using NIRCam/F150W and HST/F160W imaging (Section 4.2). Second,
we note that with HST van der Wel et al. (2014) was only able to measure the
rest-optical size evolution up to $z=2.5$, whereas with JADES we can now
extend the size measurements of quiescent galaxies beyond $z>2.5$. If we only
include $z<2.5$ galaxies in the fit, we will still find a similarly steeper
size evolution than van der Wel et al. (2014). Therefore the steeper size
evolution from our measurements is not due to including $z>2.5$ quiescent
galaxies either. The main cause for the steeper size evolution seems to be
sample selection. The analysis of van der Wel et al. (2014) included all UVJ
quiescent galaxies selected based on 3D-HST rest-frame colors, while we
removed from our analysis 36 out of 197 galaxies which do not satisfy the UVJ
criteria any more after we re-fit their SEDs with new photometry from JADES
(Section 3). If those 36 galaxies were not removed, we would have found
$\beta=-1.30\pm 0.27$ which is fully consistent with van der Wel et al.
(2014). We note that the majority – 22 out of 36 – of the removed galaxies are
at $z>1.5$, most of which are near the UVJ selection boundaries.
Interestingly, the removed galaxies generally have more extended morphologies
than those retained in our final sample, indicating a strong correlation
between the size and star formation properties of galaxies at Cosmic Noon
(Wuyts et al., 2011; Williams et al., 2010; Patel et al., 2013; van der Wel et
al., 2014; Ji & Giavalisco, 2023, just to name a few).

We now compare our results with a very recent JWST study done by van der Wel
et al. (2023) who measured the size evolution of massive quiescent galaxies
using NIRCam data from the CEERS survey. van der Wel et al. (2023) obtained
$\beta=-1.70\pm 0.10$ for $\log(M_{*}/M_{\sun})>11$, $0.5<z<2.3$ quiescent
galaxies. If we fit the size evolution only for a subsample of the same
redshift and stellar-mass ranges, we will obtain a fully consistent
$\beta\approx-1.75$. Moreover, van der Wel et al. (2023) found that the pace
of the size evolution does not depend on size proxies used (e.g., half-light
vs half-mass). This is also similar to our finding above that $\beta$ does not
depend on rest-frame wavelengths where the sizes are measured (Table 1 and
Figure 8).

Our results – and hence those of van der Wel et al. 2023 – then suggest
significant size evolution of massive quiescent galaxies beyond $z>1$, which
is in tension with some previous studies based on HST imaging: Suess et al.
(2020) and Miller et al. (2023) found that the half-mass size evolution of
quiescent galaxies is slower than their rest-optical, half-light size
evolution, suggesting that the latter is mostly due to stellar mass-to-light
radio gradients, rather than the intrinsic growth of galaxies’ sizes. In
Figure 12, for the subsample having $R_{\rm{halfmass}}$ from Suess et al.
(2020), we plot the relative difference between $R_{\rm{halfmass}}$ and
$R_{e}^{1\mu m}$ as a function of redshift. We also plot the running median
and its uncertainty estimated following the bootstrap Monte Carlo method
described earlier (Section 5.2). Up to $z\sim 1.2$, the median of the relative
difference remains roughly flat, with $R_{\rm{halfmass}}$ being larger than
$R_{e}^{1\mu m}$ by $14\pm 8\%$. Although the sample size of quiescent
galaxies having both $R_{\rm{halfmass}}$ and $R_{e}^{1\mu m}$ measurements is
limiting at $z>1.2$, we see evidence that $R_{\rm{halfmass}}$ becomes
increasingly larger than $R_{e}^{1\mu m}$ towards higher redshifts. This
explains the shallower half-mass size evolution at $z>1$ reported by those
earlier HST studies. van der Wel et al. (2023, see their Appendix A) also
attempted to address this tension regarding the half-mass size evolution at
$z>1$ by comparing the light-to-mass weighed size ratios from their
measurements with those from Suess et al. (2020). Similar to what we found in
Figure 12, van der Wel et al. 2023 showed that, while the light-to-mass
weighed size ratio from Suess et al. (2020) becomes increasingly larger at
$z>1.5$, they did not see this trend in their new measurements using CEERS
data.

We point out that the methods used to derive the half-mass sizes are
dramatically different in the literature, which can lead to inconsistent
conclusions about the size evolution of quiescent galaxies. For example, also
using HST imaging, Mosleh et al. (2017) adopted a different approach to
account for the radial gradient of the stellar mass-to-light ratio. Unlike
Suess et al. (2020) and Miller et al. (2023), Mosleh et al. (2017) found that
half-mass sizes of quiescent galaxies significantly evolve at $z>1$. In the
future, dedicated analysis of different methodologies in deriving the half-
mass sizes, and larger samples of quiescent galaxies at $z>2$ are needed to
finally address this tension between JWST and HST studies.

Figure 12: Relative difference between $R_{e}^{1\mu m}$ and the half-mass
sizes from Suess et al. (2019b) as a function of redshift. The grey dashdotted
line and shaded region show the running median the its 1$\sigma$ uncertainty.
The uncertainty of the relative difference (y-axis) is dominated by the
uncertainty of half-mass size measurements from Suess et al. (2019b).

###  5.3 Stellar mass-size Relationship of Massive Quiescent Galaxies

We now present the mass-size relationship of massive quiescent galaxies. To
characterize it, we assume the functional form of

$R_{e}^{\rm{1\mu
m}}\,[\rm{kpc}]=A\cdot\left(\frac{M_{*}}{10^{10.7}M_{\sun}}\right)^{\Gamma},$
(4)

where both the amplitude $A$ and power law index $\Gamma$ depend on redshifts
following $A=A_{0}(1+z)^{\alpha}$ and $\Gamma=\Gamma_{0}(1+z)^{\gamma}$. It
has been shown by previous studies that the mass-size relationship of
quiescent galaxies cannot be well described by a single power law where the
relationship becomes flattened below the pivot mass $\sim 10^{10.4}M_{\sun}$
(e.g. van der Wel et al., 2014; Cutler et al., 2022; van der Wel et al.,
2023). Therefore, we only fit Equation 4 for galaxies with
$M_{*}>10^{10.4}M_{\sun}$. Like before, we estimate the uncertainties of the
fitted parameters using the bootstrap Monte Carlo method: We bootstrap the
entire sample, use normal distributions to Monte Carlo resample the $M_{*}$
and $R_{e}^{\rm{1\mu m}}$ measures using their corresponding uncertainties,
and finally fit the mass-size relationship. We repeat this procedure 10000
times and use the range between 16th and 84th percentiles as $1\sigma$
uncertainties.

In Figure 13 we show the mass-size relationship. The amplitude $A$ strongly
evolves with redshift following $(1+z)^{-1.8}$, showing that at a fixed
stellar mass, quiescent galaxies are significantly larger at lower redshifts.
The slope of the mass-size relationship is $\Gamma_{0}\sim 0.5$ and it does
not strongly depend on redshift with $\gamma$ being consistent with 0, in
broad agreement with previous studies (van der Wel et al., 2014; Mowla et al.,
2019). We note, however, that the current sample size of $z>3$ quiescent
galaxies is very limiting and future larger samples are needed to check if the
slope of the mass-size relationship persists towards higher redshifts.
Finally, for completeness we also report the best-fit mass-size relationships
of the rest-frame 0.3$\micron$ and 0.5$\micron$ sizes in Table 3.

Figure 13: Stellar mass-size relationship of massive quiescent galaxies. In the left panel we show the scatter plot of stellar masses vs rest-1$\micron$ sizes, where the entire sample is divided into three groups using redshift terciles. Small dots show the individual measurements, and large squares show the median sizes of individual stellar-mass bins. For galaxies with $\log(M_{*}/M_{\sun})\geq 10.4$, we fit the redshift-dependent mass-size relationship using the functional form shown on the right (see Section 5.3 for details). Solid lines show the best-fit relationships at the median redshifts of the three groups. We also extrapolate the relation to $\log(M_{*}/M_{\sun})<10.4$, shown as dashed lines. The right panel shows the corner plots for the fitted parameters of the mass-size relationship. Like before, GS-9209 is marked with the orange square. Table 3: The best-fit mass-size relationship of massive quiescent galaxies (Equation 4). $\lambda_{\rm{rest}}$ |  $A_{0}$ [kpc] | $\alpha$ | $\Gamma_{0}$ | $\gamma$
---|---|---|---|---
$0.3\micron$ | 8.43${}_{-1.58}^{+2.10}$ | -1.78${}_{-0.30}^{+0.26}$ | 0.49${}_{-0.28}^{+0.49}$ | 0.17${}_{-0.94}^{+0.78}$
$0.5\micron$ | 5.04${}_{-0.73}^{+1.02}$ | -1.41${}_{-0.26}^{+0.22}$ | 0.48${}_{-0.20}^{+0.31}$ | 0.33${}_{-0.71}^{+0.63}$
$1\micron$ | 5.56${}_{-0.99}^{+1.27}$ | -1.81${}_{-0.30}^{+0.26}$ | 0.51${}_{-0.21}^{+0.34}$ | 0.22${}_{-0.75}^{+0.57}$

##  6 Discussion

With the deep, multi-band JWST/NIRCam imaging data from JADES, we have studied
the size evolution of a robust sample of massive quiescent galaxies with
$\log(M_{*}/M_{\sun})>10$ over $0.5<z<5$. Regardless of the rest-frame
wavelengths where the sizes are measured, we observed strong size evolution of
massive quiescent galaxies. This evolution depends on stellar mass, with more
massive quiescent galaxies having faster evolution in size. In the following,
we will discuss in detail the implications of the observed size evolution on
the formation and evolution of quiescent galaxies across cosmic time. We will
also discuss our findings with regard to the quiescent galaxies at $z>3$.

###  6.1 Physical Mechanisms Driving the Size Evolution of Quiescent Galaxies

Two broad categories of physical mechanisms can drive the strong size
evolution of quiescent galaxies without significantly increasing their stellar
masses. First, in $\Lambda$CDM cosmology, at a fixed mass galaxies formed
later are expected to be larger, the progenitor effect (Carollo et al., 2013;
Lilly & Carollo, 2016; Ji & Giavalisco, 2022). It is therefore possible that
the apparent size evolution of quiescent galaxies is driven by the addition of
newly quenched, larger galaxies at lower redshifts. Alternatively, post-
quenching growth – particularly mechanisms such as gas-poor minor mergers that
can efficiently grow sizes without significantly increasing masses – can also
explain the observed strong size evolution (e.g. Bezanson et al., 2009; Naab
et al., 2009).

Regarding the progenitor effect, there is a simple prediction for the pace of
apparent size evolution. Specifically, by construction, the Virial radius
($R_{\rm{vir}}$) of a dark matter halo is set by the threshold of density
contrast – the ratio of density to the mean density of the Universe – at the
redshift of halo collapse ($z_{\rm{form}}^{\rm{h}}$). Because the density
contrast threshold has a very weak time dependence, we then expect that, at a
fixed halo mass, $R_{\rm{vir}}$ should linearly scale with the size of the
Universe, i.e. $R_{\rm{vir}}\propto H(z_{\rm{form}}^{\rm{h}})^{-2/3}$ where
$H(z_{\rm{form}}^{\rm{h}})$ is the Hubble constant at $z_{\rm{form}}^{\rm{h}}$
(Mo et al., 2010). Given that $H(z)\sim(1+z)^{1.4}$ at Cosmic Noon $z\sim 2$,
we would then expect $R_{\rm{Vir}}\propto(1+z_{\rm{form}}^{\rm{h}})^{-1}$.
Finally, because the ratio of the galaxy $R_{e}$ to the halo $R_{\rm{Vir}}$ is
roughly constant, as supported by both theories (e.g., Mo et al., 2010) and
observations (e.g., Kravtsov, 2013; Somerville et al., 2018), the prediction
for the apparent size evolution driven by the progenitor effect is then
$R_{e}\propto(1+z_{\rm{form}}^{\rm{h}})^{-1}$.

####  6.1.1 Lower-mass quiescent galaxies with $10<\log(M_{*}/M_{\sun})<10.6$

The size evolution of $10^{10}M_{\sun}<M_{*}<10^{10.6}M_{\sun}$ quiescent
galaxies follows $(1+z)^{-1.15_{-0.22}^{+0.25}}$ (Section 5.2), which is
broadly consistent with $R_{e}\propto(1+z)^{-1}$ and hence is in line with the
progenitor effect.

Before moving forward, we clarify several things with regard to the comparison
between the observed size evolution and the prediction of the progenitor
effect. First, the prediction is made for galaxies at a fixed halo mass,
rather than a fixed stellar mass. However, we note that the stellar-mass range
considered here – $10^{10}M_{\sun}<M_{*}<10^{10.6}M_{\sun}$ – is where
galaxies’ halo masses do not change substantially over $0.5<z<5$ (see e.g.,
the left panel of Figure 7 in Behroozi et al., 2013a). If we use the
relationship between stellar mass and halo mass from Behroozi et al. (2013a,
b)888At $z\leq 5$, the data release of Behroozi et al. (2013a, b) only
contains six redshift slices, i.e. $z=0.1,1,2,3,4,5$. We therefore interpolate
those slices to derive the halo mass of a galaxy at its exact redshift value.,
we find that the halo masses of our quiescent samples at different redshifts
are very similar, with the median halo mass being $10^{12}M_{\sun}$ for the
$z>2$ subsample compared to $10^{11.9}M_{\sun}$ for the $z<2$ subsample.

Second, the halo formation theory detailed above, if applicable, suggests that
the size of galaxies should more tightly correlate with
$z_{\rm{form}}^{\rm{h}}$, the time of halo formation, than with $z$, the time
of observation. Following Ji & Giavalisco (2022), instead of using $z$, we
study the relationship between $R_{e}$ and the formation redshift,
$z_{\rm{form}}$, which is defined as the redshift when a galaxy reached its
mass-weighted age. As Figure 14 shows, the best-fit relationship is
$R_{e}^{1\mu m}\propto(1+z_{\rm{form}})^{-0.93\pm 0.27}$. Compared to
$\beta=-1.15$ obtained for the relationship between $R_{e}^{1\mu m}$ and $z$
(Section 5.2), using $z_{\rm{form}}$ makes the slope of the size evolution
closer to $-1$, i.e. the prediction of the progenitor effect. We note,
however, that using $z_{\rm{form}}$ as a proxy for $z_{\rm{form}}^{\rm{h}}$ is
purely from an empirical standpoint, and very likely there are better metrics
which however is beyond the scope of this work.

Finally, in addition to fitting the size evolution as a function of $(1+z)$,
which is widely used in the literature but arguably not the best physically
motivated form in terms of halo formation (see description above and also van
der Wel et al. 2014), we also fit the size evolution as a function of the
Hubble constant $H(z)$, i.e., $R_{e}\propto H(z)^{\beta_{H}}$. We find
$\beta_{H}=-0.87_{-0.15}^{+0.21}$ for the relationship between $R_{e}^{1\mu
m}$ and $H(z)$, and $\beta_{H}=-0.64_{-0.17}^{+0.23}$ for the relationship
between $R_{e}^{1\mu m}$ and $H(z_{\rm{form}})$. Similarly, the best-fit
$\beta_{H}$ of the relationship between $R_{e}^{1\mu m}$ and $z_{\rm{form}}$
is closer to $\beta_{H}=-2/3$ predicted by the progenitor effect.

Obviously, the interpretation of the observed slope of size evolution is
highly degenerated, as the slope might reflect either a single, dominated
physical mechanism or the interplay among several mechanisms that can mimic a
similar size evolution. Nonetheless, particularly because the slope becomes
closer to the predicted value of the progenitor effect after switching from
$z$ to $z_{\rm{form}}$, we believe that the progenitor effect may be the
primary mechanism, if not the only one, driving the apparent size evolution
observed in the quiescent galaxies with $10<\log(M_{*}/M_{\sun})<10.6$. If
this indeed is not a coincidence, our findings will then indicate that the
size of quiescent galaxies evolves at a similar pace following their dark
matter halos. Because of this, we speculate that dark matter halos play a key
role, if not the only one, in the quenching and structural evolution of
galaxies of this stellar-mass range, which has been suggested by some
theoretical works (e.g., Birnboim & Dekel, 2003; Cattaneo et al., 2006; Dekel
& Birnboim, 2006; Chen et al., 2020).

Figure 14: $R_{e}^{1\mu m}$ as a function of formation redshift
($z_{\rm{form}}$) for quiescent galaxies with $10<\log(M_{*}/M_{\sun})<10.6$.
The relationship between $R_{e}^{1\mu m}$ and $z_{\rm{form}}$ is broadly
consistent with $(1+z_{\rm{form}})^{-1}$, in line with the expectation from
the progenitor effect (Section 6.1.1).

####  6.1.2 Higher-mass quiescent galaxies with $\log(M_{*}/M_{\sun})>10.6$

Relative to the lower-mass ones, the size of higher-mass quiescent galaxies
increases faster with cosmic time (Section 5.2.1). The higher-mass quiescent
galaxies also have a steeper relationship between $R_{e}^{1\mu m}$ and
$z_{\rm{form}}$ with a slope of $-1.57\pm 0.26$. These deviate from the
prediction of the progenitor effect alone, suggesting additional mechanism(s)
contributing to the structural growth.

Higher-mass quiescent galaxies are hosted by more massive halos formed in
larger overdensities (Mo et al., 2010). Consequently, they are expected to
experience more merger events than lower-mass galaxies (e.g., Peng et al.,
2010b; Hopkins et al., 2010b; Rodriguez-Gomez et al., 2015; O’Leary et al.,
2021). Therefore, a likely explanation for their apparent steeper size
evolution is the increasing contribution of post-quenching size growth from
mergers, especially minor mergers that can grow galaxies’ sizes more
efficiently than major mergers for a fixed amount of mass to be added to the
central galaxy (e.g., Naab et al., 2009; Bezanson et al., 2009).

If post-quenching mergers do play a significant role in growing the size of
higher-mass quiescent galaxies, one would then expect that the slope of the
relationship between $R_{e}^{1\mu m}$ and $z_{\rm{form}}$ depends on galaxy
ages, as older galaxies have more time to alter their structure through
mergers. Motivated by this, we divide the higher-mass quiescent sample into
the young and old subsamples. Because galaxies at high redshifts are generally
younger, instead of making the division using the median age of the entire
sample, we first estimate the running median of stellar ages as a function of
redshift. The old (young) subsample contains galaxies older (younger) than the
median stellar age of their redshifts.

Figure 15: $R_{e}^{1\mu m}$ as a function of formation redshift
($z_{\rm{form}}$) for quiescent galaxies with $\log(M_{*}/M_{\sun})>10.6$. We
divide the sample into old (red) and young (blue) subsamples (Section 6.1.2).
The left panel shows the best-fit relationship between $R_{e}^{1\mu m}$ and
$z_{\rm{form}}$. The right panel shows the distributions of $\beta$ for the
two subsamples from our bootstrap Monte Carlo method (Section 6.1.2). The
black circle with error bars shows the median and standard deviation of
$\beta$ derived for the entire $\log(M_{*}/M_{\sun})>10.6$ sample. The
vertical, black dashed line marks $\beta=-1$, the expectation from the
progenitor effect. Old quiescent galaxies have a steeper $R_{e}^{1\mu
m}$-$z_{\rm{form}}$ relationship, while the relationship of young quiescent
galaxies is closer to $\beta=-1$, i.e. the progenitor effect.

We study the relationship between $R_{e}^{1\mu m}$ and $z_{\rm{form}}$ for the
young and old subsamples, respectively. To estimate uncertainties of the
fitted parameters, we use the bootstrap Monte Carlo method that has been
described previously (e.g., Section 5.2 and 5.3). As Figure 15 shows, although
the current sample size is still limiting, we see evidence that the
relationship between $R_{e}^{1\mu m}$ and $z_{\rm{form}}$ depends on stellar
ages. The old subsample has a steeper (more negative) slope than the young
subsample. Interestingly, although the slope ($-1.30\pm 0.56$) observed in the
young quiescent subsample is still somewhat steeper, it is broadly consistent
with $-1$ within uncertainty, and is closer to the expectation from the
progenitor effect than the old subsample. This finding is consistent with the
picture that, apart from the progenitor effect, mergers and/or post-quenching
continuous gas accretion drive additional, apparent size growth in very
massive quiescent galaxies. As they become older, massive quiescent galaxies
experience more merger events that can effectively grow their sizes by adding
mass to outskirts.

The conclusions above based on size evolution are in line with other analysis
of quiescent galaxies at Cosmic Noon. First, our conclusions are consistent
with the stacking analysis of quiescent galaxies’ surface brightness profiles
at $z\sim 2$ (see Section 3.4 of Ji & Giavalisco, 2022) and at $z\sim 0$
(Huang et al., 2013a, b), and also in broad agreement with other earlier
studies, both observationally and theoretically, where they concluded the mass
and size growths of massive quiescent galaxies have a significant contribution
from minor mergers at $z<2$ (e.g., Bezanson et al., 2009; van Dokkum et al.,
2010; Oser et al., 2012; Rodriguez-Gomez et al., 2016). Our conclusions are
further supported by a recent JWST study by Suess et al. (2023) who reported
the detection of a large population of low-mass companions within 35 kpc of
quiescent galaxies with $\log(M_{*}/M_{\sun})>10.5$ at $0.5<z<3$. Our
conclusions are also in line with recent ALMA observations for a small number
of massive quiescent galaxies at $z\sim 0.7$, where the larger-size extended
ones simultaneously show evidence for minor merging both in their SFHs and
newly accreted gas (Spilker et al., 2018; Woodrum et al., 2022).

While so far we lack the kinematics constraints for our high-redshift
quiescent galaxies, it is worth pointing out that the observed stellar-mass
dependence of the size evolution is consistent with the dynamical studies of
Early Type Galaxies (ETGs) at lower redshifts. At intermediate redshifts,
Bezanson et al. (2018) showed that ETGs at $z\sim 0.8$ become increasingly
velocity dispersion dominated as they become more massive. In addition, by
studying the excess kurtosis $h_{4}$ of the stellar velocity distribution in
$z<1$ quiescent galaxies, D’Eugenio et al. (2023) showed that $h_{4}$ – a
sensitive probe of the radial anisotropy of the stellar velocity distribution
– increases with cosmic time. At $z\sim 0$, Cappellari (2016, and references
therein) showed that the distribution of the dynamical state of ETGs is
bimodal – fast and slow rotators characterized by very different ratios of
rotational velocity to velocity dispersion. The two populations can be very
well separated using the stellar mass, with the slow rotators being more
massive than the faster rotators (see a recent review by Cappellari, 2016).
All these dynamical studies suggest that very massive galaxies after quenching
continue growing their sizes via merging with other galaxies through dynamical
friction. Because multiple incoherent merging events are expected to cause the
loss of angular momentum (e.g., Emsellem et al., 2011), very massive galaxies
eventually become verlocity dispersion dominated, i.e. slower rotators .

###  6.2 The Structure of Quiescent Galaxies at $z>3$

The deep JWST/NIRCam imaging from JADES enables us to select robust samples of
quiescent galaxies, and to measure their rest-optical and NIR morphologies
with high precision beyond $z>3$. Despite that the current sample size is very
limiting, our results suggest that the $z>3$ quiescent galaxies in our sample
have very small sizes, potentially being even smaller than extrapolating the
size evolution of lower-redshift quiescent galaxies (Section 5.2). A similar
finding is reported by Wright et al. (2023) who found evidence that $z>3$
massive quiescent galaxies in CEERS also seem to be more compact than
previously anticipated. In what follows, we take a close look at the
structural properties of the $z>3$ quiescent galaxies in this study.

####  6.2.1 Large stellar mass surface densities of the freshly quenched $z>3$
galaxies

Among eleven photometrically robust quiescent galaxies at $z>3$ in our sample,
eight have reliable $R_{e}$ estimates from single Sérsic fitting (Section
4.2). We calculate their stellar mass surface densities within the half-light
radius following $\Sigma_{e}=M_{*}/(2\pi R_{e}^{2})$.

The remaining three galaxies have best-fit $n=8$, and they have very compact
cores, as PSF spikes are clearly seen in their F444W images. While the
galaxies with $n=8$ were removed from all analysis presented above, here we
attempt to give reasonable $\Sigma_{e}$ estimates for these three galaxies as
well. Instead of using $R_{e}$ from single Sérsic fitting, we thus use the
sizes of the three galaxies from the Richardson–Lucy deconvolution approach
(Section 4.2.1).

Figure 16: Stellar mass surface density within effective radius ($\Sigma_{e}$)
vs $M_{*}$. Background contours show density distributions of UVJ-selected
star-forming (blue) and quiescent (red) galaxies at $1<z<3$ from CANDELS.
There are eleven robust quiescent galaxies at $z>3$ in our sample, among which
eight (shown as orange stars) have reliable size estimates from single Sérsic
fitting, including the galaxy GS-9209 that is highlighted with an enlarged
symbol. For the remaining three galaxies (shown as orange squares), we
calculate their $\Sigma_{e}$ using alternative size measures (see Section
6.2.1 for details). The grey horizontal dashed line marks the maximum value of
stellar mass surface density theoretically derived by Grudić et al. (2019).

We compare the $\Sigma_{e}$ of the $z>3$ quiescent galaxies with galaxies at
$1<z<3$ from CANDELS. For CANDELS galaxies, we make use of the stellar-mass
measures from Lee et al. (2018) and rest-optical sizes measured with HST/WFC3
imaging. For consistency, we thus use rest-0.5$\micron$ sizes to calculate
$\Sigma_{e}$ for the $z>3$ galaxies, although we have already checked that any
of our conclusions will not change if we use rest-1$\micron$ sizes instead. As
Figure 16 shows, the vast majority of the $z>3$ quiescent galaxies are very
compact, with $\Sigma_{e}>10^{10}\,M_{\sun}/\rm{kpc}^{2}$. Such a stellar mass
surface density is comparable to that observed in nuclear star clusters in
nearby galaxies (Norris et al., 2014), and it is also above the averaged value
found in quiescent galaxies at lower redshifts and significantly exceeds that
observed in star-forming galaxies of similar stellar masses at Cosmic Noon
(e.g. Barro et al., 2017; Ji & Giavalisco, 2023).

We note that the presence of AGN can make galaxies appear smaller than they
actually are, because the nonstellar, nucleated light from AGN can lead to a
smaller size from single Sérsic fitting than the intrinsic one of stellar
distribution (Ji et al., 2022). Despite that AGN-hosting galaxies have already
been removed from our sample through various AGN selection methods (Section
3), there could still be some missed ones because the AGN are too faint to be
identified by those selection methods. One example is GS-9209 – it shows broad
H$\alpha$ emission line but is not identified as an AGN based on SED fitting –
that we will discuss in detail in Section 6.2.2. Although we currently do not
know how many of the $z>3$ quiescent galaxies in our sample might host a
similarly faint broad-line AGN like GS-9209, we argue this does not seem to be
able to explain the strong size evolution, hence the very large $\Sigma_{e}$
of the $z>3$ quiescent galaxies for the following two main reasons. First, the
SED shapes of the quiescent galaxies all seem to be typical for
quiescent/post-starburst galaxies, with no clear evidence of hosting strong
AGN of known types. Second, the broad emission lines from AGN, if present,
likely can only affect the rest-0.5$\micron$ sizes (due to e.g., [O
III]5007Å), while they should have much less effects on the rest-1$\micron$
sizes because this wavelength range is mostly free of strong emission lines.
Yet, we still observe the very strong size evolution at rest-frame 1$\micron$.
Moreover, even for AGN hosting galaxies (e.g. X-ray or IRAC-selected AGN),
their rest-optical light is still predominantly from stars. The impact of AGN,
if not accounted for during morphological fitting, can only affect the size
measurement on the order of $\approx 10\%$ (see Section 3.2.2 and Figure 7 in
Ji et al. 2022), which is insufficient to explain the very large $\Sigma_{e}$
obtained here. Therefore, we argue that the very large $\Sigma_{e}$ found here
are most likely due to the very compact stellar-mass distributions in the
$z>3$ quiescent galaxies.

The strong correlation between the mass surface density and star-formation
properties found at low redshifts has long been suggested as an indication of
the causal link between the buildup of central regions of galaxies and
quenching (Kauffmann et al., 2003; Franx et al., 2008; Cheung et al., 2012;
Barro et al., 2017; Lee et al., 2018; Ji & Giavalisco, 2023, just to name a
few). However, because the most quiescent galaxies observed at $z<3$ are
already old, with a typical stellar age of $>1$ Gyr (Carnall et al., 2019;
Tacchella et al., 2022; Ji & Giavalisco, 2022), such a timescale is comparable
to/longer than some post-quenching processes required to alter galaxy
structures (e.g., see Figure 10 in Tacconi et al., 2020). This makes it
difficult to definitively answer whether the large mass surface densities
found in quiescent galaxies are the cause or the consequence of quenching, or
a by-product of an as-yet unidentified physical process. Compared to the
lower-redshift ones, the $z>3$ quiescent galaxies presented here are much
younger, with stellar ages ranging from 100$-$700 Myr and a typical value of
500 Myr. These freshly quenched galaxies are already compact, suggesting the
possible coeval formation of dense central regions and quenching at $z>3$
which is consistent with a recent JWST study of $3<z<4.5$ galaxies with
forming stellar cores (Ji et al., 2023).

####  6.2.2 GS-9209 at $z_{spec}=4.658$

The most remarkable $z>3$ quiescent galaxy in our sample is GS-9209, one of
the earliest massive quiescent galaxy known so far having spectroscopically
confirmed $z_{\rm{spec}}=4.658$ with JWST/NIRSpec observations (Carnall et
al., 2023b). The NIRSpec spectrum of GS-9209 is very similar to typical post-
starburst galaxies, showing prominent stellar absorption features of A-type
stars. Despite that AGN is not needed to fit GS-9209’s broad-band SED, a broad
H$\alpha$ emission line with FWHM $\sim 10^{4}$ km/s is observed in the
NIRSpec spectrum, indicating that the galaxy actually hosts a faint AGN
(Carnall et al., 2023b). Moreover, as Figure 16 shows, the galaxy has an
extreme $\Sigma_{e}$, namely that it is comparable to the maximum allowed
stellar mass surface density set by stellar feedback (Grudić et al., 2019).

With deep NIRCam images from JADES, we measure the multi-band morphologies for
GS-9209. The results are shown in Figure 17. Carnall et al. (2023b) reported a
size of $215\pm 20$ pc using the NIRCam F210M (rest-0.4$\micron$) image from
JEMS (Williams et al., 2023), which is consistent with our measurement where
we obtain the half-light semi-major radius of $223\pm 16$ pc (or $186\pm 13$
pc for circularized size $R_{e}$) in F200W. As the top panel of Figure 17
shows, overall the size of GS-9209 decreases with wavelength – its light
distribution becomes more compact at longer wavelength. The $R_{e}$ of the
galaxy is $160\pm 11$ pc in F277W (rest-0.5$\micron$) and $168\pm 12$ pc in
F444W (rest-1$\micron$), which are $2\sim 3$ times smaller than the sizes
predicted by using our best-fit size evolution (Section 5.2). In addition,
shown in the middle panel of Figure 17, the Sérsic index $n$ also changes with
wavelength. The morphology of GS-9209 gets more spheroidal at longer
wavelengths. The best-fit Sérsic index is $n=1.8\pm 0.4$ in F090W
(rest-1500Å), whereas it becomes $n=3.3\pm 0.4$ in F444W (rest-1$\micron$).
These immediately show the presence of radial color gradient in GS-9209. The
galaxy is more extended in rest-UV and optical than in rest-NIR, either due to
the radial gradient of dust attenuation where the central regions of the
galaxy either have more dust, or older stellar populations, or both.

Finally, we also note the sudden change of the light distribution of GS-9209
in F356W where the H$\alpha$ emission of the galaxy is covered. The galaxy is
more nucleated (Sérsic $n=5.2\pm 0.7$) in F356W than in adjacent filters, as
the AGN H$\alpha$ emission comes from the very central parts of the galaxy. In
the future, we may be able to recover some faint high-redshift AGN candidates
missed by photometric selection methods through morphological analysis of deep
NIRCam images, once the galaxies’ redshifts are well constrained.

Figure 17: Multi-band morphologies of GS-9209 in NIRCam images from JADES. The
black, vertical dashed line marks H$\alpha$ emission, and the shaded grey
region corresponds to FWHM $=10,000$ km/s measured by Carnall et al. (2023b).
Generally speaking, from rest-UV to NIR, GS-9209 becomes more compact (smaller
$R_{e}$) and spheroidal (larger Sérsic $n$). Evidence of broad H$\alpha$
affecting the light distribution in F356W is also observed.

##  7 Summary

In this work, we studied the size evolution of massive quiescent galaxies with
$\log(M_{*}/M_{\sun})>10$ at $0.5<z<5$. With multi-band deep NIRCam imaging
from the JADES survey in the GOODS-South field, we characterized the size
evolution of massive quiescent galaxies from UV through NIR at three different
wavelengths, namely rest-frame 0.3, 0.5 and 1$\micron$. We found that

  1. 1.

Regardless of wavelengths, the size of massive quiescent galaxies strongly
evolves with cosmic time following $R_{e}\propto(1+z)^{-1.3}$. The pace of
size evolution depends on stellar mass, with the evolution being faster for
more massive galaxies. Specifically, the size evolution follows
$R_{e}\propto(1+z)^{-1.1}$ for lower-mass quiescent galaxies with
$10<\log(M_{*}/M_{\sun})<10.6$ , while it follows $R_{e}\propto(1+z)^{-1.7}$
for higher-mass quiescent galaxies with $\log(M_{*}/M_{\sun})>10.6$ (Section
5.2.1).

  2. 2.

As a result of negative stellar mass-to-light ratio, the size of quiescent
galaxies decreases with wavelength. Relative to their rest-1$\micron$ sizes,
on average quiescent galaxies are larger by 45% at rest-0.3$\micron$ and 15%
at rest-0.5$\micron$, showing that the stellar mass distribution of quiescent
galaxies is more compact than the rest-UV and optical light distributions.
Evidence that younger quiescent galaxies have shallower color gradients than
older ones is also seen, but requires future larger quiescent samples to test
(Section 5.2.2).

  3. 3.

At a fixed stellar mass, the size of quiescent galaxies is significantly
larger at lower redshifts. The slope of the relationship between stellar mass
and size only very weakly, if at all, depends on redshift, although we stress
that the current sample size of $z>3$ quiescent galaxies is very limiting and
future larger samples are needed to check if the slope remains unchanged
towards higher redshifts (Section 5.3).

We discussed possible physical mechanisms responsible for the observed size
evolution. In particular, we studied the relationship between $R_{e}$ and
$z_{\rm{form}}$, the formation redshift of galaxies. We showed that

  1. 1.

For lower-mass ($10<\log(M_{*}/M_{\sun})<10.6$) quiescent galaxies, the slope
of the relationship between $R_{e}$ and $z_{\rm{form}}$ is consistent with
$-1$, the prediction of the progenitor effect. Thus we think that the
progenitor effect may play a significant role in driving the apparent size
evolution of quiescent galaxies of this stellar-mass range (Section 6.1.1).

  2. 2.

For higher-mass ($\log(M_{*}/M_{\sun})>10.6$) quiescent galaxies, the
relationship between $R_{e}$ and $z_{\rm{form}}$ depends on stellar age. For
young quiescent galaxies, the slope of the $R_{e}$-$z_{\rm{form}}$
relationship is consistent with $-1$ within uncertainty (though somewhat
steeper), whereas for old quiescent galaxies the relationship is much steeper
than the prediction of the progenitor effect alone. These are in line with the
picture that, in addition to the progenitor effect, mergers and/or post-
quenching continuous gas accretion drive additional, apparent size growth in
very massive quiescent galaxies (Section 6.1.2).

Notwithstanding the very limited sample size, we studied the structure of
massive quiescent galaxies at $z>3$ in this work. We saw evidence that these
quiescent galaxies are even smaller than anticipated using the best-fit size
evolution derived for lower-redshift quiescent galaxies (Section 5.2.1).
Moreover, we showed that these freshly quenched quiescent galaxies, with a
typical stellar age of 500 Myr compared to $>1$ Gyr observed in quiescent
galaxies at lower redshifts, are already very compact, with $\Sigma_{e}\gtrsim
10^{10}M_{\sun}/\rm{kpc}^{2}$, suggesting the possible coeval formation of
dense central regions and quenching at $z>3$ (Section 6.2). We also studied in
detail the multi-band morphology for GS-9209 (Carnall et al., 2023b), one of
the earliest massive quiescent galaxies known so far with
$z_{\rm{spec}}=4.658$. From rest-UV to NIR, we found that GS-9209 becomes
increasingly smaller and more compact, and its light profile becomes more
spheroidal, showing that the radial color gradient is already present in this
earliest massive quiescent galaxy (Section 6.2.2).

Finally, as a closing remark, we note that the current high-redshift quiescent
sample, e.g. at $z>3$, is still small. Explicitly and better constraining the
size evolution at $z>3$ will await larger samples selected from JWST imaging,
and this motivates future analysis based on robustly JWST-selected high-
redshift samples (e.g., Alberts et al., 2023) over larger sky areas. We also
stress that characterizing the size evolution really is just the first step to
constrain the formation pathways of the first quiescent galaxies in the
Universe. The structure of the earliest quiescent galaxies is likely much more
complex. Some hints have already been seen from our Figure 5: Sizes from the
Richardson–Lucy deconvolution approach are larger than those from the single
Sérsic fitting, suggesting multi-component structures are already present in
quiescent galaxies at $z>3$. With the superb NIR imaging from JWST, we will be
able to study such complex structural properties in the earliest quiescent
galaxies in unprecedented detail.

ZJ, BDJ, BER, ZC, CD, DJE, EE, KH, JMH, GR, MR, FS and CNAW acknowledge
support from JWST/NIRCam contract to the University of Arizona, NAS5-02015.
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. BER acknowledges
support from JWST Program 3215. SA acknowledges support from the JWST Mid-
Infrared Instrument (MIRI) Science Team Lead, grant 80NSSC18K0555, from NASA
Goddard Space Flight Center to the University of Arizona. WB, FDE, TJL, RM, LS
and JW acknowledge support by the Science and Technology Facilities Council
(STFC), by the ERC through Advanced Grant 695671 “QUENCH”. KB acknowledges
support in part by the Australian Research Council Centre of Excellence for
All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number
CE170100013. AJB and JC acknowledge funding from the “FirstGalaxies” Advanced
Grant from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (Grant agreement No. 789056).
SC acknowledges support by European Union’s HE ERC Starting Grant No.
101040227 - WINGS. ECL acknowledges support of an STFC Webb Fellowship
(ST/W001438/1). DJE is supported as a Simons Investigator. RH acknowledges
that funding for this research was provided by the Johns Hopkins University,
Institute for Data Intensive Engineering and Science (IDIES). RM acknowledges
support by the UKRI Frontier Research grant RISEandFALL, and funding from a
research professorship from the Royal Society. HÜ gratefully acknowledges
support by the Isaac Newton Trust and by the Kavli Foundation through a
Newton-Kavli Junior Fellowship. This work made use of the lux supercomputer at
UC Santa Cruz which is funded by NSF MRI grant AST 1828315, as well as the
High Performance Computing (HPC) resources at the University of Arizona which
is funded by the Office of Research Discovery and Innovation (ORDI), Chief
Information Officer (CIO), and University Information Technology Services
(UITS).

##  Appendix A SED fitting with different SFHs

We test our SED fitting with different assumptions of SFHs. Specifically, in
Figure 18, we compare the physical quantities derived using our fiducial SED
assumptions, i.e. nonparametric SFHs with the Dirichlet prior (Section 4.1),
with those derived using nonparametric SFHs with the continuity prior. These
two nonparametric SFH priors are very similar, except that the continuity
prior disfavor sudden changes of SFRs in adjacent lookback time bins (see Leja
et al. 2019 for details). Similarly, in Figure 19, we compare the our fiducial
SED fitting results with those derived using parametric delayed-tau SFHs.
Stellar masses, redshifts and rest-frame UVJ colors agree very well among
different SFH assumptions. Systematic offsets exist for the SFR measures for
quiescent galaxies, as expected for SED fitting with photometric data only
(Conroy, 2013).

Figure 18: Comparison of physical properties derived from nonparametric SFHs
with the Dirichlet prior ($x$-axis) and nonparametric SFHs with the continuity
prior ($y$-axis). Figure 19: Similar to Figure 18, but now we compare the
physical properties from nonparametric Dirichlet SFHs with those from
parametric delayed-tau SFHs.

## References

  * Alberts et al. (2023) Alberts, S., Williams, C. C., Helton, J. M., et al. 2023, arXiv e-prints, arXiv:2312.12207.  https://arxiv.org/abs/2312.12207
  * Anderson & King (2000) Anderson, J., & King, I. R. 2000, PASP, 112, 1360, doi: 10.1086/316632
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Barro et al. (2017) Barro, G., Faber, S. M., Koo, D. C., et al. 2017, ApJ, 840, 47, doi: 10.3847/1538-4357/aa6b05
  * Behroozi et al. (2013a) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013a, ApJ, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Behroozi et al. (2013b) —. 2013b, ApJ, 762, L31, doi: 10.1088/2041-8205/762/2/L31
  * Bell & de Jong (2001) Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212, doi: 10.1086/319728
  * Belli et al. (2019) Belli, S., Newman, A. B., & Ellis, R. S. 2019, ApJ, 874, 17, doi: 10.3847/1538-4357/ab07af
  * Bezanson et al. (2009) Bezanson, R., van Dokkum, P. G., Tal, T., et al. 2009, ApJ, 697, 1290, doi: 10.1088/0004-637X/697/2/1290
  * Bezanson et al. (2018) Bezanson, R., van der Wel, A., Pacifici, C., et al. 2018, ApJ, 858, 60, doi: 10.3847/1538-4357/aabc55
  * Birnboim & Dekel (2003) Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349, doi: 10.1046/j.1365-8711.2003.06955.x
  * Bradley et al. (2023) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2023, astropy/photutils: 1.8.0, 1.8.0, Zenodo, doi: 10.5281/zenodo.7946442
  * Brammer et al. (2012) Brammer, G. B., van Dokkum, P. G., Franx, M., et al. 2012, ApJS, 200, 13, doi: 10.1088/0067-0049/200/2/13
  * Byler et al. (2017) Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44, doi: 10.3847/1538-4357/aa6c66
  * Cappellari (2016) Cappellari, M. 2016, ARA&A, 54, 597, doi: 10.1146/annurev-astro-082214-122432
  * Carnall et al. (2019) Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2019, MNRAS, 490, 417, doi: 10.1093/mnras/stz2544
  * Carnall et al. (2020) Carnall, A. C., Walker, S., McLure, R. J., et al. 2020, MNRAS, 496, 695, doi: 10.1093/mnras/staa1535
  * Carnall et al. (2023a) Carnall, A. C., McLeod, D. J., McLure, R. J., et al. 2023a, MNRAS, 520, 3974, doi: 10.1093/mnras/stad369
  * Carnall et al. (2023b) Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2023b, arXiv e-prints, arXiv:2301.11413, doi: 10.48550/arXiv.2301.11413
  * Carollo et al. (2013) Carollo, C. M., Bschorr, T. J., Renzini, A., et al. 2013, ApJ, 773, 112, doi: 10.1088/0004-637X/773/2/112
  * Cassata et al. (2013) Cassata, P., Giavalisco, M., Williams, C. C., et al. 2013, ApJ, 775, 106, doi: 10.1088/0004-637X/775/2/106
  * Cattaneo et al. (2006) Cattaneo, A., Dekel, A., Devriendt, J., Guiderdoni, B., & Blaizot, J. 2006, MNRAS, 370, 1651, doi: 10.1111/j.1365-2966.2006.10608.x
  * Charlot & Fall (2000) Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718, doi: 10.1086/309250
  * Chen et al. (2020) Chen, Z., Faber, S. M., Koo, D. C., et al. 2020, ApJ, 897, 102, doi: 10.3847/1538-4357/ab9633
  * Cheung et al. (2012) Cheung, E., Faber, S. M., Koo, D. C., et al. 2012, ApJ, 760, 131, doi: 10.1088/0004-637X/760/2/131
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393, doi: 10.1146/annurev-astro-082812-141017
  * Conroy & Gunn (2010) Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833, doi: 10.1088/0004-637X/712/2/833
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486
  * Cutler et al. (2022) Cutler, S. E., Whitaker, K. E., Mowla, L. A., et al. 2022, ApJ, 925, 34, doi: 10.3847/1538-4357/ac341c
  * Daddi et al. (2005) Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680, doi: 10.1086/430104
  * Damjanov et al. (2011) Damjanov, I., Abraham, R. G., Glazebrook, K., et al. 2011, ApJ, 739, L44, doi: 10.1088/2041-8205/739/2/L44
  * Davari et al. (2014) Davari, R., Ho, L. C., Peng, C. Y., & Huang, S. 2014, ApJ, 787, 69, doi: 10.1088/0004-637X/787/1/69
  * Dekel & Birnboim (2006) Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2, doi: 10.1111/j.1365-2966.2006.10145.x
  * D’Eugenio et al. (2023) D’Eugenio, F., van der Wel, A., Piotrowska, J. M., et al. 2023, MNRAS, 525, 2789, doi: 10.1093/mnras/stad800
  * Donley et al. (2012) Donley, J. L., Koekemoer, A. M., Brusa, M., et al. 2012, ApJ, 748, 142, doi: 10.1088/0004-637X/748/2/142
  * Donnari et al. (2019) Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485, 4817, doi: 10.1093/mnras/stz712
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8, doi: 10.3847/0067-0049/222/1/8
  * Eisenstein et al. (2023) Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv e-prints, arXiv:2306.02465, doi: 10.48550/arXiv.2306.02465
  * Emsellem et al. (2011) Emsellem, E., Cappellari, M., Krajnović, D., et al. 2011, MNRAS, 414, 888, doi: 10.1111/j.1365-2966.2011.18496.x
  * Falcón-Barroso et al. (2011) Falcón-Barroso, J., Sánchez-Blázquez, P., Vazdekis, A., et al. 2011, A&A, 532, A95, doi: 10.1051/0004-6361/201116842
  * Fan et al. (2010) Fan, L., Lapi, A., Bressan, A., et al. 2010, ApJ, 718, 1460, doi: 10.1088/0004-637X/718/2/1460
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Franx et al. (2008) Franx, M., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2008, ApJ, 688, 770, doi: 10.1086/592431
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001, doi: 10.1088/1538-3873/acd1b5
  * Giavalisco et al. (2004) Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600, L93, doi: 10.1086/379232
  * Glazebrook et al. (2023) Glazebrook, K., Nanayakkara, T., Schreiber, C., et al. 2023, arXiv e-prints, arXiv:2308.05606, doi: 10.48550/arXiv.2308.05606
  * Grudić et al. (2019) Grudić, M. Y., Hopkins, P. F., Quataert, E., & Murray, N. 2019, MNRAS, 483, 5548, doi: 10.1093/mnras/sty3386
  * Guo et al. (2013) Guo, Y., Ferguson, H. C., Giavalisco, M., et al. 2013, ApJS, 207, 24, doi: 10.1088/0067-0049/207/2/24
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hopkins et al. (2009) Hopkins, P. F., Bundy, K., Murray, N., et al. 2009, MNRAS, 398, 898, doi: 10.1111/j.1365-2966.2009.15062.x
  * Hopkins et al. (2010a) Hopkins, P. F., Murray, N., Quataert, E., & Thompson, T. A. 2010a, MNRAS, 401, L19, doi: 10.1111/j.1745-3933.2009.00777.x
  * Hopkins et al. (2010b) Hopkins, P. F., Bundy, K., Croton, D., et al. 2010b, ApJ, 715, 202, doi: 10.1088/0004-637X/715/1/202
  * Huang et al. (2013a) Huang, S., Ho, L. C., Peng, C. Y., Li, Z.-Y., & Barth, A. J. 2013a, ApJ, 766, 47, doi: 10.1088/0004-637X/766/1/47
  * Huang et al. (2013b) —. 2013b, ApJ, 768, L28, doi: 10.1088/2041-8205/768/2/L28
  * Ilbert et al. (2013) Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55, doi: 10.1051/0004-6361/201321100
  * Illingworth et al. (2016) Illingworth, G., Magee, D., Bouwens, R., et al. 2016, arXiv e-prints, arXiv:1606.00841, doi: 10.48550/arXiv.1606.00841
  * Ji & Giavalisco (2022) Ji, Z., & Giavalisco, M. 2022, ApJ, 935, 120, doi: 10.3847/1538-4357/ac7f43
  * Ji & Giavalisco (2023) —. 2023, ApJ, 943, 54, doi: 10.3847/1538-4357/aca807
  * Ji et al. (2022) Ji, Z., Giavalisco, M., Kirkpatrick, A., et al. 2022, ApJ, 925, 74, doi: 10.3847/1538-4357/ac3837
  * Ji et al. (2018) Ji, Z., Giavalisco, M., Williams, C. C., et al. 2018, ApJ, 862, 135, doi: 10.3847/1538-4357/aacc2c
  * Ji et al. (2020) Ji, Z., Giavalisco, M., Vanzella, E., et al. 2020, ApJ, 888, 109, doi: 10.3847/1538-4357/ab5fdc
  * Ji et al. (2023) Ji, Z., Williams, C. C., Tacchella, S., et al. 2023, arXiv e-prints, arXiv:2305.18518, doi: 10.48550/arXiv.2305.18518
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22, doi: 10.3847/1538-4365/abef67
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 54, doi: 10.1046/j.1365-8711.2003.06292.x
  * Kravtsov (2013) Kravtsov, A. V. 2013, ApJ, 764, L31, doi: 10.1088/2041-8205/764/2/L31
  * Kubo et al. (2018) Kubo, M., Tanaka, M., Yabe, K., et al. 2018, ApJ, 867, 1, doi: 10.3847/1538-4357/aae3e8
  * Lee et al. (2018) Lee, B., Giavalisco, M., Whitaker, K., et al. 2018, ApJ, 853, 131, doi: 10.3847/1538-4357/aaa40f
  * Leja et al. (2019) Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3, doi: 10.3847/1538-4357/ab133c
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170, doi: 10.3847/1538-4357/aa5ffe
  * Leja et al. (2022) Leja, J., Speagle, J. S., Ting, Y.-S., et al. 2022, ApJ, 936, 165, doi: 10.3847/1538-4357/ac887d
  * Lilly & Carollo (2016) Lilly, S. J., & Carollo, C. M. 2016, ApJ, 833, 1, doi: 10.3847/0004-637X/833/1/1
  * Lower et al. (2020) Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33, doi: 10.3847/1538-4357/abbfa7
  * Luo et al. (2017) Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2, doi: 10.3847/1538-4365/228/1/2
  * Lustig et al. (2021) Lustig, P., Strazzullo, V., D’Eugenio, C., et al. 2021, MNRAS, 501, 2659, doi: 10.1093/mnras/staa3766
  * Lyu et al. (2022) Lyu, J., Alberts, S., Rieke, G. H., & Rujopakarn, W. 2022, ApJ, 941, 191, doi: 10.3847/1538-4357/ac9e5d
  * Madau (1995) Madau, P. 1995, ApJ, 441, 18, doi: 10.1086/175332
  * Maltby et al. (2018) Maltby, D. T., Almaini, O., Wild, V., et al. 2018, MNRAS, 480, 381, doi: 10.1093/mnras/sty1794
  * McGaugh & Schombert (2014) McGaugh, S. S., & Schombert, J. M. 2014, AJ, 148, 77, doi: 10.1088/0004-6256/148/5/77
  * Merlin et al. (2021) Merlin, E., Castellano, M., Santini, P., et al. 2021, A&A, 649, A22, doi: 10.1051/0004-6361/202140310
  * Miller et al. (2023) Miller, T. B., van Dokkum, P., & Mowla, L. 2023, ApJ, 945, 155, doi: 10.3847/1538-4357/acbc74
  * Miller et al. (2022) Miller, T. B., Whitaker, K. E., Nelson, E. J., et al. 2022, ApJ, 941, L37, doi: 10.3847/2041-8213/aca675
  * Mo et al. (2010) Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution
  * Mosleh et al. (2020) Mosleh, M., Hosseinnejad, S., Hosseini-ShahiSavandi, S. Z., & Tacchella, S. 2020, ApJ, 905, 170, doi: 10.3847/1538-4357/abc7cc
  * Mosleh et al. (2017) Mosleh, M., Tacchella, S., Renzini, A., et al. 2017, ApJ, 837, 2, doi: 10.3847/1538-4357/aa5f14
  * Mowla et al. (2019) Mowla, L., van der Wel, A., van Dokkum, P., & Miller, T. B. 2019, ApJ, 872, L13, doi: 10.3847/2041-8213/ab0379
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18, doi: 10.1088/0004-637X/777/1/18
  * Naab et al. (2009) Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJ, 699, L178, doi: 10.1088/0004-637X/699/2/L178
  * Nanayakkara et al. (2022) Nanayakkara, T., Glazebrook, K., Jacobs, C., et al. 2022, arXiv e-prints, arXiv:2212.11638, doi: 10.48550/arXiv.2212.11638
  * Nelson et al. (2016) Nelson, E. J., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2016, ApJ, 828, 27, doi: 10.3847/0004-637X/828/1/27
  * Newman et al. (2012) Newman, A. B., Ellis, R. S., Bundy, K., & Treu, T. 2012, ApJ, 746, 162, doi: 10.1088/0004-637X/746/2/162
  * Norris et al. (2014) Norris, M. A., Kannappan, S. J., Forbes, D. A., et al. 2014, MNRAS, 443, 1151, doi: 10.1093/mnras/stu1186
  * O’Leary et al. (2021) O’Leary, J. A., Moster, B. P., Naab, T., & Somerville, R. S. 2021, MNRAS, 501, 3215, doi: 10.1093/mnras/staa3746
  * Oser et al. (2012) Oser, L., Naab, T., Ostriker, J. P., & Johansson, P. H. 2012, ApJ, 744, 63, doi: 10.1088/0004-637X/744/1/63
  * Pacifici et al. (2023) Pacifici, C., Iyer, K. G., Mobasher, B., et al. 2023, ApJ, 944, 141, doi: 10.3847/1538-4357/acacff
  * Papovich et al. (2018) Papovich, C., Kawinwanichakij, L., Quadri, R. F., et al. 2018, ApJ, 854, 30, doi: 10.3847/1538-4357/aaa766
  * Patel et al. (2013) Patel, S. G., van Dokkum, P. G., Franx, M., et al. 2013, ApJ, 766, 15, doi: 10.1088/0004-637X/766/1/15
  * Peebles (1993) Peebles, P. J. E. 1993, Principles of Physical Cosmology, doi: 10.1515/9780691206721
  * Peng et al. (2010a) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010a, AJ, 139, 2097, doi: 10.1088/0004-6256/139/6/2097
  * Peng et al. (2010b) Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010b, ApJ, 721, 193, doi: 10.1088/0004-637X/721/1/193
  * Perrin et al. (2014a) Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014a, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, ed. J. Oschmann, Jacobus M., M. Clampin, G. G. Fazio, & H. A. MacEwen, 91433X, doi: 10.1117/12.2056689
  * Perrin et al. (2014b) Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014b, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, ed. J. Oschmann, Jacobus M., M. Clampin, G. G. Fazio, & H. A. MacEwen, 91433X, doi: 10.1117/12.2056689
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Prato et al. (2012) Prato, M., Cavicchioli, R., Zanni, L., Boccacci, P., & Bertero, M. 2012, A&A, 539, A133, doi: 10.1051/0004-6361/201118681
  * Richardson (1972) Richardson, W. H. 1972, Journal of the Optical Society of America (1917-1983), 62, 55
  * Rieke & the JADES Collaboration (2023) Rieke, M., & the JADES Collaboration. 2023, arXiv e-prints, arXiv:2306.02466, doi: 10.48550/arXiv.2306.02466
  * Robotham et al. (2017) Robotham, A. S. G., Taranu, D. S., Tobar, R., Moffett, A., & Driver, S. P. 2017, MNRAS, 466, 1513, doi: 10.1093/mnras/stw3039
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 2015, MNRAS, 449, 49, doi: 10.1093/mnras/stv264
  * Rodriguez-Gomez et al. (2016) Rodriguez-Gomez, V., Pillepich, A., Sales, L. V., et al. 2016, MNRAS, 458, 2371, doi: 10.1093/mnras/stw456
  * Schreiber et al. (2015) Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74, doi: 10.1051/0004-6361/201425017
  * Shen et al. (2003) Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978, doi: 10.1046/j.1365-8711.2003.06740.x
  * Skelton et al. (2014) Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24, doi: 10.1088/0067-0049/214/2/24
  * Somerville et al. (2018) Somerville, R. S., Behroozi, P., Pandya, V., et al. 2018, MNRAS, 473, 2714, doi: 10.1093/mnras/stx2040
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15, doi: 10.1088/0067-0049/214/2/15
  * Spilker et al. (2018) Spilker, J., Bezanson, R., Barišić, I., et al. 2018, ApJ, 860, 103, doi: 10.3847/1538-4357/aac438
  * Straatman et al. (2015) Straatman, C. M. S., Labbé, I., Spitler, L. R., et al. 2015, ApJ, 808, L29, doi: 10.1088/2041-8205/808/1/L29
  * Suess et al. (2019a) Suess, K. A., Kriek, M., Price, S. H., & Barro, G. 2019a, ApJ, 885, L22, doi: 10.3847/2041-8213/ab4db3
  * Suess et al. (2019b) —. 2019b, ApJ, 877, 103, doi: 10.3847/1538-4357/ab1bda
  * Suess et al. (2020) —. 2020, ApJ, 899, L26, doi: 10.3847/2041-8213/abacc9
  * Suess et al. (2022) Suess, K. A., Bezanson, R., Nelson, E. J., et al. 2022, ApJ, 937, L33, doi: 10.3847/2041-8213/ac8e06
  * Suess et al. (2023) Suess, K. A., Williams, C. C., Robertson, B., et al. 2023, arXiv e-prints, arXiv:2307.14209, doi: 10.48550/arXiv.2307.14209
  * Tacchella et al. (2015) Tacchella, S., Carollo, C. M., Renzini, A., et al. 2015, Science, 348, 314, doi: 10.1126/science.1261094
  * Tacchella et al. (2022) Tacchella, S., Conroy, C., Faber, S. M., et al. 2022, ApJ, 926, 134, doi: 10.3847/1538-4357/ac449b
  * Tacconi et al. (2020) Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157, doi: 10.1146/annurev-astro-082812-141034
  * Toft et al. (2007) Toft, S., van Dokkum, P., Franx, M., et al. 2007, ApJ, 671, 285, doi: 10.1086/521810
  * Tortora et al. (2010) Tortora, C., Napolitano, N. R., Cardone, V. F., et al. 2010, MNRAS, 407, 144, doi: 10.1111/j.1365-2966.2010.16938.x
  * Trujillo et al. (2007) Trujillo, I., Conselice, C. J., Bundy, K., et al. 2007, MNRAS, 382, 109, doi: 10.1111/j.1365-2966.2007.12388.x
  * Valentino et al. (2023) Valentino, F., Brammer, G., Gould, K. M. L., et al. 2023, arXiv e-prints, arXiv:2302.10936, doi: 10.48550/arXiv.2302.10936
  * van der Wel et al. (2011) van der Wel, A., Rix, H.-W., Wuyts, S., et al. 2011, ApJ, 730, 38, doi: 10.1088/0004-637X/730/1/38
  * van der Wel et al. (2012) van der Wel, A., Bell, E. F., Häussler, B., et al. 2012, ApJS, 203, 24, doi: 10.1088/0067-0049/203/2/24
  * van der Wel et al. (2014) van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28, doi: 10.1088/0004-637X/788/1/28
  * van der Wel et al. (2023) van der Wel, A., Martorano, M., Haussler, B., et al. 2023, arXiv e-prints, arXiv:2307.03264, doi: 10.48550/arXiv.2307.03264
  * van Dokkum et al. (2008) van Dokkum, P. G., Franx, M., Kriek, M., et al. 2008, ApJ, 677, L5, doi: 10.1086/587874
  * van Dokkum et al. (2010) van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018, doi: 10.1088/0004-637X/709/2/1018
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Whitaker et al. (2013) Whitaker, K. E., van Dokkum, P. G., Brammer, G., et al. 2013, ApJ, 770, L39, doi: 10.1088/2041-8205/770/2/L39
  * Whitaker et al. (2019) Whitaker, K. E., Ashas, M., Illingworth, G., et al. 2019, ApJS, 244, 16, doi: 10.3847/1538-4365/ab3853
  * Williams et al. (2017) Williams, C. C., Giavalisco, M., Bezanson, R., et al. 2017, ApJ, 838, 94, doi: 10.3847/1538-4357/aa662f
  * Williams et al. (2023) Williams, C. C., Tacchella, S., Maseda, M. V., et al. 2023, arXiv e-prints, arXiv:2301.09780, doi: 10.48550/arXiv.2301.09780
  * Williams et al. (2009) Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879, doi: 10.1088/0004-637X/691/2/1879
  * Williams et al. (2010) Williams, R. J., Quadri, R. F., Franx, M., et al. 2010, ApJ, 713, 738, doi: 10.1088/0004-637X/713/2/738
  * Woodrum et al. (2022) Woodrum, C., Williams, C. C., Rieke, M., et al. 2022, ApJ, 940, 39, doi: 10.3847/1538-4357/ac9af7
  * Wright et al. (2023) Wright, L., Whitaker, K. E., Weaver, J. R., et al. 2023, arXiv e-prints, arXiv:2311.05394, doi: 10.48550/arXiv.2311.05394
  * Wuyts et al. (2011) Wuyts, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011, ApJ, 742, 96, doi: 10.1088/0004-637X/742/2/96
  * Xue et al. (2010) Xue, Y. Q., Brandt, W. N., Luo, B., et al. 2010, ApJ, 720, 368, doi: 10.1088/0004-637X/720/1/368