# Size - Stellar Mass Relation and Morphology of Quiescent Galaxies at $z\geq
3$ in Public $JWST$ Fields

Kei Ito JSPS Research Fellow (PD<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Astronomy, School of Science, The
University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan
Francesco Valentino European Southern Observatory, Karl-Schwarzschild-Str. 2,
D-85748 Garching bei Munchen, Germany  Cosmic Dawn Center (DAWN), Denmark
Gabriel Brammer  Cosmic Dawn Center (DAWN), Denmark  Niels Bohr Institute,
University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark  Andreas
L. Faisst Caltech/IPAC, 1200 E. California Blvd. Pasadena, CA 91125, USA
Steven Gillman Cosmic Dawn Center (DAWN), Denmark  DTU-Space, Technical
University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark  Carlos
Gómez-Guijarro Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM,
91191, Gif-sur-Yvette, France  Katriona M. L. Gould  Cosmic Dawn Center
(DAWN), Denmark  Niels Bohr Institute, University of Copenhagen, Jagtvej 128,
DK-2200 Copenhagen N, Denmark  Kasper E. Heintz  Cosmic Dawn Center (DAWN),
Denmark  Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200
Copenhagen N, Denmark  Olivier Ilbert  Aix Marseille Univ, CNRS, CNES, LAM,
Marseille, France  Christian Kragh Jespersen Department of Astrophysical
Sciences, Princeton University, Princeton, NJ 08544, USA  Vasily Kokorev
Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700AV
Groningen, The Netherlands  Cosmic Dawn Center (DAWN), Denmark  Mariko Kubo
Astronomical Institute, Tohoku University, Aoba-ku, Sendai 980-8578, Japan
Georgios E. Magdis  Cosmic Dawn Center (DAWN), Denmark  DTU-Space, Technical
University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark  Niels
Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N,
Denmark  Conor McPartland  Cosmic Dawn Center (DAWN), Denmark  Niels Bohr
Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N,
Denmark  Masato Onodera  Department of Astronomical Science, The Graduate
University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo,
181-8588, Japan  Subaru Telescope, National Astronomical Observatory of Japan,
National Institutes of Natural Sciences (NINS), 650 North A’ohoku Place, Hilo,
HI 96720, USA  Francesca Rizzo  Cosmic Dawn Center (DAWN), Denmark  Niels Bohr
Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N,
Denmark  Masayuki Tanaka  National Astronomical Observatory of Japan, 2-21-1
Osawa, Mitaka, Tokyo, 181-8588, Japan  Department of Astronomical Science, The
Graduate University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka,
Tokyo, 181-8588, Japan  Sune Toft  Cosmic Dawn Center (DAWN), Denmark  Niels
Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N,
Denmark  Aswin P. Vijayan  Cosmic Dawn Center (DAWN), Denmark  DTU-Space,
Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark
John R. Weaver Department of Astronomy, University of Massachusetts, Amherst,
MA 01003, USA  Katherine E. Whitaker  Department of Astronomy, University of
Massachusetts, Amherst, MA 01003, USA  Cosmic Dawn Center (DAWN), Denmark
Lillian Wright Department of Physics and Astronomy, Seoul National University,
1 Gwanak-ro, Gwanak-gu, Seoul, South Korea  Department of Astronomy,
University of Massachusetts, Amherst, MA 01003, USA

###### Abstract

We present the results of a systematic study of the rest-frame optical
morphology of quiescent galaxies at $z\geq 3$ using the Near-Infrared Camera
(NIRCam) onboard $JWST$. Based on a sample selected by $UVJ$ color or $NUVUVJ$
color, we focus on 26 quiescent galaxies with
$9.8<\log{(M_{\star}/M_{\odot})}<11.4$ at $2.8<z_{\rm phot}<4.6$ with publicly
available $JWST$ data. Their sizes are constrained by fitting the Sérsic
profile to all available NIRCam images. We see a negative correlation between
the observed wavelength and the size in our sample and derive their size at
the rest-frame $0.5\,{\rm\mu m}$ taking into account this trend. Our quiescent
galaxies show a significant correlation between the rest-frame $0.5\,{\rm\mu
m}$ size and the stellar mass at $z\geq 3$. The analytical fit for them at
$\log{(M_{\star}/M_{\odot})}>10.3$ implies that our size - stellar mass
relations are below those at lower redshifts, with the amplitude of $\sim
0.6\,{\rm kpc}$ at $M_{\star}=5\times 10^{10}\,M_{\odot}$. This value agrees
with the extrapolation from the size evolution of quiescent galaxies at $z<3$
in the literature, implying that the size of quiescent galaxies increases
monotonically from $z\sim 3-5$. Our sample is mainly composed of galaxies with
bulge-like structures according to their median Sérsic index and axis ratio of
$n\sim 3-4$ and $q\sim 0.6-0.8$, respectively. On the other hand, there is a
trend of increasing fraction of galaxies with low Sérsic index, suggesting
$3<z<5$ might be the epoch of onset of morphological transformation with a
fraction of very notable disky quenched galaxies.

Galaxy evolution (594); Galaxy quenching (2040); High-redshift galaxies (734);
Galaxy radii (617); Quenched galaxies (2016)

††facilities: JWST††software: Astropy (Robitaille et al., 2013; Price-Whelan
et al., 2018), emcee (Foreman-Mackey et al., 2013), galfit (Peng et al.,
2010), lmfit (Newville et al., 2014), Matplotlib (Hunter, 2007), numba (Lam et
al., 2015), numpy (Harris et al., 2020), pandas (McKinney, 2010), photutils
(Bradley et al., 2022), SciPy (Virtanen et al., 2020), WebbPSF (Perrin et al.,
2012, 2014)

##  1 Introduction

One of the most fundamental ways to characterize galaxies is to look at how
they look like. The size and morphology correlate with the formation history
of galaxies, including major or minor mergers (e.g., Bezanson et al., 2009),
and their host halos. Thus, tracing them across cosmic time provides us with
insights into the evolution of galaxies (e.g., van der Wel et al., 2014, and
references therein).

The relationship between the galaxy half-light size in the rest-frame optical
and stellar mass (hereafter called size-mass relation) has been examined up to
$z\sim 3$ (e.g., Shen et al., 2003; van der Wel et al., 2014; Faisst et al.,
2017; Mowla et al., 2019b; Kawinwanichakij et al., 2021; Cutler et al., 2022),
showing that the size is larger for higher stellar mass galaxies. The slope of
the size-mass relation is steeper for quiescent galaxies than that for star-
forming galaxies, which could be due to the dry minor mergers (see the
discussion in e.g., van Dokkum et al., 2015). Moreover, these size-mass
relations suggest that the size of quiescent galaxies drastically evolves as a
function of redshift. For example, van der Wel et al. (2014) show that the
size of quiescent galaxies at fixed stellar mass evolves as
$\propto(1+z)^{-1.48}$ at $z<3$, whereas that of star-forming galaxies evolves
as $\propto(1+z)^{-0.75}$. Various phenomena can cause this redshift evolution
of size. Firstly, newly quenched galaxies at lower redshift having a larger
size similar to star-forming galaxies can cause an apparent size evolution,
which is called progenitor bias (e.g., Cassata et al., 2013; Carollo et al.,
2013; Barro et al., 2013). The dry minor mergers also effectively increase the
size of quiescent galaxies (e.g., Oser et al., 2012; Nipoti et al., 2012). In
addition, this trend might be only in the half-light size. Suess et al. (2019)
shows that the half-mass size does not show as strong redshift evolution
compared to the half-light size mentioned above (see also Mosleh et al.,
2017). This implies that the observed light profile is different dependent on
the wavelength. To minimize its effect, some studies derive the half-light
radius in the same rest-frame wavelength in comparing the different redshifts
by using images spanning a range of observed-frame wavelengths (e.g.,
Kawinwanichakij et al., 2021; Nedkova et al., 2021).

Some quiescent galaxies at $z\sim 2-3$ are reported as disk galaxies or spiral
galaxies (e.g., Toft et al., 2017; Newman et al., 2018; Fudamoto et al.,
2022). In contrast, the majority of quiescent galaxies have a high Sérsic
index and axis ratio even at high redshift, which implies that they have
bulge-like shapes (e.g., Straatman et al., 2015; Marsan et al., 2019; Lustig
et al., 2021). Thus, there could be growing morphological variation for
quiescent galaxies at high redshift.

Recently, massive quiescent galaxies have been found even at $z\gtrsim 3-4$
(Glazebrook et al., 2017; Schreiber et al., 2018; Tanaka et al., 2019; Forrest
et al., 2020; Valentino et al., 2020; D’Eugenio et al., 2021; Nanayakkara et
al., 2022; Carnall et al., 2023a, b). The morphology of some of these galaxies
are also explored. For example, Kubo et al. (2018) stack the $K$-band images
assisted by Adaptive Optics (AO) and derive the average rest-frame optical
size of quiescent galaxies at $z\sim 4$, which includes a spectroscopically
confirmed quiescent galaxy at $z=4.01$ reported in Tanaka et al. (2019) and
Valentino et al. (2020). Also, the Hubble Space Telescope ($HST$)/WFC3 images
constrain the morphological properties of some $z\geq 3$ spectroscopically
confirmed quiescent galaxies (Gobat et al., 2012; Saracco et al., 2020; Lustig
et al., 2021; Esdaile et al., 2021).

However, a statistical discussion of the structural properties did not be
carried out for quiescent galaxies at $z\geq 3$, including the constraint on
their size-mass relation and its redshift evolution to $z\sim 0$. There are
four underlying reasons. One is the requirement for high-resolution images to
resolve these galaxies at $\geq 2\,{\rm\mu m}$; the rest-frame wavelength of
$0.5\,{\rm\mu m}$ at $z\geq 3$, often used in the size-mass relation at lower
redshift, is not observable with $HST$. For $z>3$ galaxies, $HST$ can only
probe the wavelength of $<0.4\,{\rm\mu m}$ in the rest frame. Other
approaches, such as using AO, lead to expensive observation and a severely
limited and incomplete sample. Secondly, since the method for selecting
quiescent galaxies differs among the aforementioned literature (e.g., based on
the specific star formation rate or colors, also see Table 4 in Valentino et
al., 2023, for the summary of the literature), comparisons between different
samples are not straightforward. Thirdly, assembling a large sample of this
rare galaxy population requires large and deep imaging, especially in near-
infrared. Finally, constraining the wavelength dependence of the size of
individual objects had been difficult. If there is a correlation between the
structural properties and the wavelength for quiescent galaxies at $z\geq 3$,
as at low redshift, the measurement with a single identical band could cause a
systematic bias due to the different redshift.

The James Webb Space Telescope ($JWST$) is poised to revolutionize this topic.
$JWST$’s high-resolution imaging at near-infrared, particularly at $\lambda>2\
\mu$m, allows us to probe the emission from the stellar populations of
galaxies at $z\geq 3$ and accurately measure their morphological properties in
the rest-frame optical. Moreover, utilizing its multi-band information across
a wide range of wavelengths, we can evaluate the wavelength dependence of the
sizes, as done at lower redshift. Indeed, several studies derive structural
properties of distant quiescent galaxies with the Near-Infrared Camera
(NIRCam) on $JWST$, such as for a quiescent galaxy at $z=4.658$ (Carnall et
al., 2023b), and a statistical sample of galaxies at $1\leq z\leq 2.5$,
including quiescent ones (Suess et al., 2022). The size evolution is also
explored for star-forming galaxies at higher redshift (e.g., Ono et al., 2022;
Yang et al., 2022; Gómez-Guijarro et al., 2023).

In this paper, we systematically investigate the structural properties of a
sample of 26 quiescent galaxies at $z\geq 3$ from Valentino et al. (2023)
(hereafter denoted as V23) by using the $JWST$/NIRCam images. V23 select
quiescent galaxies in fields with publicly available $JWST$ data based on
$U-V$ and $V-J$ color (Williams et al., 2009) or $NUV-U$, $U-V$, and $V-J$
color (Gould et al., 2023). This paper derives the rest-frame optical size-
stellar mass relation and quantifies their typical morphological properties,
focusing on the Sérsic index and axis ratio in particular. This paper is
organized as follows. The sample selection for this study is summarized in
Section 2. In Section 3, we summarize the method of morphological fitting and
deriving the rest-optical size. Section 4 reports the size-mass relations of
our quiescent galaxies and the analytical fit to them. The morphological
properties, i.e., Sérsic index and axis ratio, are discussed in Section 5.
Lastly, we discuss and summarize our results in Section 6. We assume the
$\Lambda$CDM cosmology with $H_{0}=70\ {\rm km\ s^{-1}Mpc^{-1}}$,
$\Omega_{m}=0.3$, and $\Omega_{\Lambda}=0.7$. The magnitude is based on the AB
magnitude system (Oke & Gunn, 1983). We assume the initial mass function of
Chabrier (2003).

##  2 Sample Selection

This study uses the sample of quiescent galaxies at $z\geq 3$ selected from
$JWST$ data, constructed in V23. Here we give a brief overview of the sample
selection and refer the reader to V23 for more details.

V23 process the all publicly available $JWST$ imaging of fields spanning $\sim
145\,{\rm arcmin^{2}}$ taken during the first three months and extract sources
therein using Source Extractor Python (sep) (Barbary, 2016). Their photometric
redshift, rest-frame colors, and stellar mass are estimated by eazy-py
(Brammer et al., 2008; Brammer, 2021) with $JWST$ and $HST$ photometry. We
confirm their photometric redshift agrees well with the spectroscopic redshift
determined from the archival observations (G. Brammer et al., in preparation).
In the CEERS field, we achieve a $\sigma_{NMAD}=0.0268$ in $\delta z/(1+z)$
for the spectroscopic sample. We select galaxies at $z\geq 3$ from their
photometric redshift.

Two color selections for quiescent galaxies are employed. One is based on the
traditional $UVJ$ diagram. V23 select galaxies located in the box of quiescent
galaxies defined in Williams et al. (2009). We here do not consider the
uncertainty of the color. This selection corresponds to the “standard”
selection argued in V23. The other is based on the $NUV-U$, $U-V$, and $V-J$
color diagram, hereafter called the $NUVUVJ$ diagram. This is recently
proposed by Gould et al. (2023), which constructs a Gaussian Mixture Model for
the quiescent galaxy selection at high redshift (Gould, 2023). Gould et al.
(2023) shows that this selection not only can reduce the contamination
compared to the $UVJ$ diagram but also can select recently quenched galaxies
more efficiently, which can be abundant at $z>3$. Based on the $NUVUVJ$
diagram, this model returns the “probability of being quiescent”, $P_{\rm Q}$.
This study selects quiescent galaxies based on the $P_{\rm Q,50\%}>0.7$, where
$P_{\rm Q,50\%}$ is the median of $P_{Q}$ derived from bootstrapping. We note
that these selections of the $UVJ$ and $NUVUVJ$ colors are the most
conservative ones not to be affected by the contaminants due to the
uncertainty in color. These $UVJ$ and $NUVUVJ$ color selections give 34 and 18
quiescent galaxies, respectively.

We further impose three selections to accurately derive morphological
properties. Firstly, we require galaxies to be detected in more than three
bands of $JWST$ with the signal-to-noise ratio of $S/N>50$ in $0.5\arcsec$
diameter aperture photometry. This criterion is motivated by a simulation
using mock Sérsic profiles summarized in Section 3.4 and ensures that we have
reliable and precise structural parameter measurements (i.e., effective
radius, Sérsic index) in multiple bands. As described below, loosening this
$S/N$ criterion will result in including more systematically biased output
values. We here retain 30 and 18 galaxies selected by $UVJ$ or $NUVUVJ$,
respectively. Second, we remove X-ray detections to avoid the contribution
from active galactic nuclei (AGNs) to their images, which results in removing
a galaxy (#2718 in V23 ID) detected in the AEGIS-X Deep survey (Nandra et al.,
2015) and a galaxy (#2552 in V23 ID) detected in the Chandra observation of
Tzanavaris et al. (2014). They are both selected only by $NUVUVJ$, and this
leads to 30 and 16 galaxies selected by $UVJ$ or $NUVUVJ$, respectively. We
note that none of the objects in our sample were listed in existing radio
catalogs. Lastly, we further do not use two quiescent galaxies from the
SMACS0723 field to avoid overestimating the sizes due to the strong
gravitational lensing.

In total, we use 28 quiescent galaxies selected by at least either $UVJ$
selection or $NUVUVJ$ selection. All of these galaxies are selected by $UVJ$,
and 14 galaxies are selected by $NUVUVJ$. They are at $2.83\leq z\leq 4.63$
with $9.77\leq\log{(M_{\star}/M_{\odot})}\leq 11.39$ according to eazy-py.
Figure 1 summarizes the RGB images of the long wavelength of NIRCam for these
galaxies. Their location in the $UVJ$ and $NUVUVJ$ diagram are shown in Figure
2. We refer the reader to the supplement data of V23 for their SED (Valentino,
2023). The sample includes galaxies in CEERS (ERS #1345, PI: S. Finkelstein,
Bagley et al., 2023), Stephan’s Quintet (ERO #2736, Pontoppidan et al., 2022),
PRIMER (PRIMER, GO #1837), North Ecliptic Pole Time-Domain Field (GTO, #2738,
PI: R. Windhorst), J1235 (COM/NIRCam #1063, PI: B.Sunnquiest), and SGAS1723
(as a part of TEMPLATES, ERS #1355, PI: J. Rigby).

Figure 1: Color images of all quiescent galaxies used in this study. F277W,
F356W, and F444W images are used for blue, green, and red, respectively. All
images are $4.0\arcsec\times 4.0\arcsec$ centered on each quiescent galaxy.
They are sorted according to their redshift. At the bottom of each figure, the
ID of galaxies, which is identical to those of V23, and the method which
selects it as quiescent one are shown. We can see that most of them have
compact profiles. Figure 2: Left panel: $UVJ$ diagram (Williams et al., 2009).
Squares represent quiescent galaxies selected by the $UVJ$ only, whereas
circles represent those selected by both $UVJ$ and $NUVUVJ$. They are colored
according to their photometric redshift, and the size of markers is larger for
more massive galaxies. The gray background represents the distribution of
galaxies with $\log{(M_{\star}/M_{\odot})}>9.5$ at $3<z<5$. Right panel:
$NUVUVJ$ diagram (Gould et al., 2023). The meanings of the markers are the
same as in the left panel.

##  3 Morphological Fitting

To obtain the structural parameters, we fit an analytical function to the
$JWST$/NIRCam images. Though there could be asymmetrical features or two
components, such as bulge and disk, we simply fit the single Sérsic profile
(Sérsic, 1963) following the literature (e.g., van der Wel et al., 2014; Mowla
et al., 2019b; Nedkova et al., 2021, and others) to fairly compare them. We
note that the reduced chi-square of the fitting described later implies that
this assumption is reasonable.

###  3.1 Input Images

We use all available images of $JWST$/NIRCam from the F115W filter to the
F444W filter (see Table 3 of V23 for the full list of available filters on
each field). We use the cutouts of $4.0\arcsec\times 4.0\arcsec$ of each
filter as input object images. The variance images generated by the Grizli
pipeline (Brammer & Matharu, 2021; Brammer et al., 2022) are used to make
sigma images. The sigma image can sometimes be overestimated or
underestimated, thus they are renormalized by a correction factor which makes
the sigma image of the images without any sources equal to its fluctuation.
Object masks are generated using the photutils module. Objects more than
$0.4\arcsec$ away from the target objects are masked out.

###  3.2 Point Spread Function

The point spread function (PSF) is an essential factor in conducting the
morphological fitting precisely. Previous studies use the theoretical PSFs
(e.g., Suess et al., 2022) based on the WebbPSF software (Perrin et al., 2012,
2014), those based on the software with the real images (e.g., Zhuang & Shen,
2023), such as PSFEx (Bertin, 2011), or natural stars (e.g., Ding et al.,
2022; Yang et al., 2022). This study utilizes the theoretical PSFs from
WebbPSF smoothed to match their radial flux profile of natural stars for each
field. This approach copes with the possible difference in the profile between
the theoretical PSFs and the natural stars reported recently (Ding et al.,
2022; Ono et al., 2022) and is free from the noises which exist in the stacked
star images.

The theoretical PSFs are first generated for each filter with WebbPSF. Here,
in-flight sensing measurements are used by retrieving the Optical Path
Difference files at the nearest dates to the observation date. The obtained
PSFs are rotated to match the position angle. Next, we derive the radial flux
profile of natural stars. Stars are selected from the source catalog of each
field by focusing on the unresolved objects based on the FLUX_RADIUS and with
the small ellipticity from sep, which is a similar approach to that employed
in PSFEx. The detailed selection is summarized in Appendix A. We then smooth
the PSFs from WebbPSF with the Moffat profile to match its one-dimensional
radial profile with that of the median-stacked natural stars. The full width
at half maximum of the used PSF is $0.06\arcsec-0.15\arcsec$, which is larger
at a longer wavelength. Figure 12 shows the one-dimensional radial profiles of
smoothed PSFs used in this study in CEERS as an example.

###  3.3 Fitting

In fitting the Sérsic profile to images, we use Galfit (Peng et al., 2002,
2010). It uses real object images, sigma images, mask images, and PSFs as
inputs and returns the best-fit parameters. We first fit the F277W image of
each galaxy, which shows the light profile at around the rest-frame
$0.5\,{\rm\mu m}$ similar to van der Wel et al. (2014) for our galaxies
($\lambda_{\rm rest}=0.49-0.72\,{\rm\mu m}$ depending on their redshift). In
the fitting, the magnitude $m$, effective radius $r_{\rm eff}$, Sérsic index
$n$, axis ratio $q$, position angle PA, and the central pixel coordinates are
set as the free parameters. We allow these parameters in the following range
to avoid the unrealistic cases: $-3\,{\rm mag}<m-m_{\rm SEP}<3\,{\rm mag}$,
$0.1\,{\rm pixel}<r_{\rm eff}<500\,{\rm pixel}$, $0.2<n<8$, and $0.0001<q<1$,
where $m_{\rm SEP}$ is the Kron magnitude from the sep. In fitting the images
in the other filters, we fix the axis ratio and position angle as the best-fit
value in the F277W images as in Nedkova et al. (2021) since they are not
likely to depend on the wavelength. This enables us to derive the effective
radius and the Sérsic index in a consistent manner in the different
wavelengths. We get the Sérsic index in the edge of the parameter range
($n=8$) for some galaxies. In that case, we refit the images of all filters
assuming $n=4$. The sky background of the images is also set to be a free
parameter, and surrounding objects within $0.4\arcsec$ from the target
galaxies are simultaneously fit by the Sérsic profile. Sources that are more
than $0.4\arcsec$ away are masked out.

Figure 3 shows an example of the result of the fitting for all available
images of a galaxy, which is the typical one of our sample with the redshift
and the stellar mass close to the median values. Appendix B also summarizes
the fitting results in the F200W filter and F277W filter and the best-fit
parameters in the F277W filter of all galaxies used in Section 4 and 5 (see
Section 3.5 for the selection). Some galaxies have compact residuals with both
positive and negative pixels. They are often seen in $HST$ images of compact
galaxies with high Sérsic index, too (e.g., Marsan et al., 2019; Lustig et
al., 2021; Forrest et al., 2022). Most of the residuals are just 10% of the
flux of the original images at most. The median reduced chi-square of the
fitting for targets is $\chi^{2}_{\nu}=1.0-1.2$ in each filter, suggesting
that the assumption of the single Sérsic profile is overall good.

To estimate realistic uncertainties of the obtained parameters, the PSF-
convolved Sérsic profile of the best-fit parameters for each galaxy is first
inserted into a noise image for each filter and each field. The noise images
are created by randomly cutting out the real images without any sources. We
apply Galfit to this image in the same manner as for the real image and obtain
the best fit. This procedure is repeated for 100 different noise images, and
the uncertainty of each parameter is taken as the range between the 16th and
84th values of this distribution. We find that the original uncertainty
reported by Galfit is 2-5 times more underestimated depending on wavelength
than that derived by this procedure.

Figure 3: A Galfit fit of CEERS #7995 as an example, which is $UVJ$-selected
quiescent galaxy at $z=3.21$ with $\log{(M_{\star}/M_{\odot})}=10.54$ in CEERS
field. The observed object image, mask image, model image, and residual image
are shown from left to right panels in each row. The longer wavelength images
are shown toward the bottom row. All images are $4.0\arcsec\times 4.0\arcsec$.
The images are scaled by the noise ($\pm 15\sigma$) of object images in each
filter.

###  3.4 Bias and quality in the fitting

The quality of the Galfit best-fit model depends on the signal-to-noise ratio
of the images (e.g., Figure 6 in van der Wel et al., 2012). Thus, the
different brightness and different profile of galaxies will have a different
quality in the fitting. Here, we evaluate the bias and the quality of the
fitting using mock galaxy images following the Sérsic profile.

Firstly, 3000 different single Sérsic profiles with various ranges of
parameters are prepared. Parameters for each profile are randomly selected in
the following range: $21\,{\rm mag}<m<29\,{\rm mag}$, $0.002\arcsec<r_{\rm
eff}<2\arcsec$, and $0.1<q<1$. The position angle is set to be random. Since
we are focusing on quiescent galaxies, which have bulge-like structures
supported by the literature (e.g., Lustig et al., 2021) and Section 5.1 of
this paper, we assume Sérsic index as $n=3-5$. Next, after convolving with the
corresponding PSF, we insert them into 30 different noise images for each. In
this section, the CEERS F277W image is used as an example. We note that the
results do not significantly change if we use images in different filters. In
total, we generate 90000 mock Sérsic profile images. We fit the suite of mock
galaxies using Galfit in the same manner as for the real galaxies and compare
the obtained parameters with the input parameters.

Figure 4 shows the relative uncertainty of the effective radius, Sérsic index,
and axis ratio. As an example, we here focus on the trend at $\sim 26.3$ mag,
which corresponds to the $S/N=50$ detection in the CEERS F277W image (see V23
for the depth of images). Firstly, the effective radii are slightly
overestimated in smaller input sizes (20% in $0.002\arcsec<r_{\rm
eff}<0.02\arcsec$) and slightly underestimated in larger input sizes (20% in
$0.2\arcsec<r_{\rm eff}<2\arcsec$). The Sérsic index is underestimated
($10-40\%$, equal to $\delta n\sim 0.4-1.6$ for $n=3-5$) regardless of the
input size and axis ratio. The axis ratio is well reproduced for
$0.02\arcsec<r_{\rm eff}<0.2\arcsec$, but underestimated in larger or smaller
sizes, in particular larger axis ratio ($q>0.4$). Fainter (brighter) objects
are expected to have more (less) significant scatter and bias in all
parameters.

As described in the following sections, our sample has $r_{\rm eff}\sim 0.6$
kpc ($\sim 0.1\arcsec$) and $q\sim 0.6-0.8$. Therefore, this test implies that
our sample with our signal-to-noise ratio cut $S/N>50$ is typically expected
to have only negligible offset in the effective radii and axis ratio ($<10\%$)
and a slight offset in the Sérsic index ($<30\%$).

Figure 4: Median relative uncertainties of the estimated parameters (effective
radius: left, Sérsic index: middle, axis ratio: right) derived through mock
galaxy images as a function of the magnitude. The top, middle, and bottom
panel corresponds to the case with the input axis ratio of $0.1<q<0.4$,
$0.4<q<0.7$, and $0.7<q<1.0$, respectively. Different colors correspond to
different input effective radii, i.e., $0.002\arcsec<r_{\rm eff}<0.02\arcsec$,
$0.02\arcsec<r_{\rm eff}<0.2\arcsec$, and $0.2\arcsec<r_{\rm eff}<2\arcsec$.
The shaded regions correspond to the ranges between the 16th and 84th values.
The vertical line in each panel corresponds to the $S/N=50$ in the CEERS F277W
image.

###  3.5 Wavelength dependence of the size, rest-frame optical sizes

It is known that the size of galaxies has a dependency on the observed
wavelength at $z<3$ (e.g., van der Wel et al., 2014; Kawinwanichakij et al.,
2021; Suess et al., 2022). Our sample spans a wide redshift range; thus,
obtaining the effective radius in the same rest-frame wavelength for all
galaxies is important to make a fair comparison with the literature at lower
redshift.

We obtain the effective radius in the rest-frame $0.5\,{\rm\mu m}$ and the
wavelength dependence of the effective radius by fitting the measured
effective radii in the different wavelengths to the following form:

$\log{r_{\rm eff}(\lambda_{\rm rest})}=\gamma\log{\left(\frac{\lambda_{\rm
rest}}{0.5\,{\rm\mu m}}\right)}+\log{r_{\rm eff,0.5\mu m}},$ (1)

where $\gamma$ and $r_{\rm eff,0.5\mu m}$ correspond to the strength of the
wavelength dependence of the effective radius ($\gamma=\Delta\log{r_{\rm
eff}}/\Delta\log{\lambda_{\rm rest}}$) and the effective radius at the rest-
frame $0.5\,{\rm\mu m}$. $\lambda_{\rm rest}$ is the wavelength in the rest
frame according to the photometric redshift. We note that the uncertainty in
photometric redshift (see Section 2) induces only negligible impact on the
rest-frame wavelength $\lambda_{\rm rest}$ ($\sim 3\%$). Following the
literature, we use the effective radius along the semi-major axis as the proxy
of the galaxy size. In the fitting, we use the effective radii of the images
in all bands where the target source is detected with $S/N>50$. This detection
threshold is motivated by the mock simulation described in Section 3.4. In
addition, the best effective radii at the edge of the parameter range (i.e.,
$r_{\rm eff}=0.1$ pixel or $500$ pixel) are not used for deriving the
effective radius at the rest-frame $0.5\,{\rm\mu m}$ and the wavelength
dependence of the size. We here require that at least three bands are used.
The number of data points in the fitting depends on the field and the target’s
brightness, which spans from three to seven.

We successfully determine $\gamma$ and $r_{\rm eff,0.5\mu m}$ of 26 quiescent
galaxies, which we use in the following analysis. Two galaxies (#7912 and
#362) are removed here since they are too compact, whose best effective radii
are 0.1 pixel ($0.004\arcsec$) in all wavelengths at least longer than $2.7\
\mu$m. Among 26 galaxies, all of them are selected by the $UVJ$ diagram, and
13 of them are selected by the $NUVUVJ$ diagram. We remind the reader that
Appendix B summarizes the fitting result and the best-fit parameters of all of
them.

Figure 5 shows the strength of the wavelength dependence of the effective
radius $\gamma$ as a function of the stellar mass. Most quiescent galaxies
have a negative correlation between the wavelength and the size, i.e., smaller
sizes in longer wavelengths. The median values of $UVJ$-selected and
$NUVUVJ$-selected sample is $\gamma=-0.38_{-0.08}^{+0.08}$ and
$\gamma=-0.33_{-0.11}^{+0.10}$, respectively. This negative correlation
between wavelength and the size agrees with the lower redshift results (e.g.,
van der Wel et al., 2014; Kawinwanichakij et al., 2021; Suess et al., 2022).
van der Wel et al. (2014) report the average strength of the wavelength
dependence of the effective radius of quiescent galaxies at $0<z<2$ as
$\gamma=-0.25$, less significant than our value. This could be due to the
redshift evolution or the different stellar mass distribution between our
sample and van der Wel et al. (2014). Given the small size of our sample, we
do not attempt to model the correlation between the stellar mass and the
strength of the wavelength dependence of the effective radius as reported in
Kawinwanichakij et al. (2021), but we note that most of our sample is in good
agreement with the relationship at the highest redshift bin ($0.8<z<1.0$) in
Kawinwanichakij et al. (2021) with some outliers with the stronger negative
gradient at lower redshift and lower mass.

Figure 5: Wavelength dependence of the effective radius
$\gamma=\Delta\log{r_{\rm eff}}/\Delta\log{\lambda_{\rm rest}}$ as a function
of stellar mass for $UVJ$ (left) and $NUVUVJ$ (right) selected quiescent
galaxies. Each marker is colored according to its redshift. Large squares show
the median of each sample. A dashed line is the relation of quiescent galaxies
at $0.8<z<1.0$ in Kawinwanichakij et al. (2021).

##  4 Rest-frame Optical Size-mass relation

###  4.1 Overall Trends

We show quiescent galaxies at $z\geq 3$ in the size-mass plane in Figure 6.
The effective radii at the rest-frame $0.5\,{\rm\mu m}$ obtained in Section
3.5 are used. First of all, we can see that many quiescent galaxies of our
sample have effective radii of less than $1$ kpc, which is significantly
smaller than those of star-forming galaxies at $z<3$ (e.g., van der Wel et
al., 2014). Their effective radii appear to correlate with their stellar mass,
i.e., more massive galaxies have larger radii. This is consistent with the
size-mass relation of quiescent galaxies as in the lower redshift studies
(e.g., van der Wel et al., 2014; Mowla et al., 2019b), which will be discussed
in the following section.

We see that several galaxies deviate from the trend. Some of the less massive
galaxies ($\log{(M_{\star}/M_{\odot})}<10.5$) have $r_{\rm eff}>2\,{\rm kpc}$
in the $UVJ$-selected case (#788, #7399, #8967). This may mean that these
galaxies are in a period of morphological transformation from star-forming
galaxies larger than quiescent galaxies or the contaminants of star-forming
galaxies. On the other hand, we can not rule out that this is due to the
difficulty in fitting the Sérsic profile to these galaxies. For example, #788
has an elongated structure, which can be a sign of a merger (Figure 13). Also,
#7399 can be affected by nearby sources. We note that removing these galaxies
from the sample does not affect our conclusions. In contrast, there are
galaxies with significantly small sizes ($r_{\rm eff}<0.1\,{\rm kpc}$, #185
and #1427 in PRIMER). AGNs can cause the apparent compact morphology, but we
do not have clear evidence that this is due to their contribution because they
are not detected in X-ray (X-UDS, Kocevski et al., 2018). The limiting flux of
X-UDS survey is $1.4\times 10^{-16}\,{\rm erg/cm^{2}/s}$ in the $0.5-2\,{\rm
keV}$. This corresponds to the limiting luminosity in the rest-frame
$0.5-10\,{\rm keV}$ of $L_{X}=(1.5-1.6)\times 10^{43}\,{\rm erg/s}$ at their
redshift, which implies that they are at least not likely powerful unobscured
X-ray AGNs.

We note that our stellar mass is robust against the codes we use. The stellar
mass of our sample measured with eazy-py is compared to that measured with
FAST++111https://github.com/cschreib/fastpp. They are in good agreement with
the offset of
$\log{(M_{\rm\star,EAZY}/M_{\odot})}-\log{(M_{\rm\star,FAST++}/M_{\odot})}=0.069_{-0.066}^{+0.075}\
(0.067_{-0.027}^{+0.033})$ for our $UVJ$ ($NUVUVJ$)-selected quiescent
galaxies, respectively. van der Wel et al. (2014), Mowla et al. (2019b) and
Nedkova et al. (2021), whose size-mass relation will be compared with ours,
use the stellar mass measured with FAST. Thus, this supports that there is no
significant bias in comparing our size-mass relation with them.

Figure 6: Size-mass relation of $UVJ$-selected (left) and $NUVUVJ$-selected
(right) quiescent galaxies at $z\geq 3$. The effective radii are measured at
the rest-frame $0.5\,\mu$m. Circles represent quiescent galaxies from this
study, color-coded by their photometric redshift. The red line and shaded
region are the best fit of the size-mass relation of quiescent galaxies from
this study and its $1\sigma$ uncertainty, respectively. The data points
located in the gray hatched region ($\log{(M_{\star}/M_{\odot})}<10.3$) are
not used for the analytical fit. The solid and dashed black line is the best
fit of the size-mass relation of early-type (i.e., quiescent) and late-type
(i.e., star-forming) galaxies at $2.5<z<3.0$ in van der Wel et al. (2014),
respectively.

###  4.2 Analytic fits

Motivated by the appearance of the existence of size-mass relation, we attempt
to analytically fit our sample following the literature (van der Wel et al.,
2014; Mowla et al., 2019b; Nedkova et al., 2021; Kawinwanichakij et al.,
2021). Firstly, following van der Wel et al. (2014), the effective radii
distribution is assumed to be a log-normal distribution $N(\log{R_{\rm
eff}},\sigma_{\log{R_{\rm eff}}})$, where $\log{R_{\rm eff}}$ is the mean of
the effective radius and $\sigma_{\log{R_{\rm eff}}}$ is the intrinsic scatter
of the distribution. We assume that the mean of the effective radius $R_{\rm
eff}$ depends on the stellar mass as:

$R_{\rm eff}/{\rm kpc}=A(\frac{M_{\star}}{5\times
10^{10}M_{\odot}})^{\alpha},$ (2)

where $A$ is the effective radius at $M_{\star}=5\times 10^{10}\,M_{\odot}$
and $\alpha$ is the slope of the size-mass relation. We consider the
uncertainty $\delta r_{\rm eff}$ of the observed effective radius $r_{\rm
eff}$ by assuming it follows a Gaussian distribution. The $\delta r_{\rm eff}$
includes not only the observed uncertainty of $r_{\rm eff}$, but also the
uncertainty of the stellar mass, assuming $0.15$ dex, which is the same value
assumed in van der Wel et al. (2014) and roughly comparable to the scatter of
eazy-py (see the discussion in Section 4.1 and Gould et al., 2023). The latter
is included by converting it to the uncertainty of the effective radius by
assuming the typical slope ($\alpha=0.6$) of the size-mass relation at lower
redshift (e.g., Mowla et al., 2019b).

Following van der Wel et al. (2014), we calculate the total likelihood for
three parameters ($A$, $\alpha$, and $\sigma_{\log{R_{\rm eff}}}$) as

$\mathcal{L}=\Sigma\ln\left[W\cdot\langle N(\log{R_{\rm
eff}},\sigma_{\log{R_{\rm eff}}}),N(\log{r_{\rm eff}},\delta\log{r_{\rm
eff}})\rangle\right],$ (3)

where $W$ is the weight to avoid the dominance of lower-mass galaxies. We use
the inverse of the stellar mass function of quiescent galaxies at $3<z<3.5$ in
Weaver et al. (2022) as the weight. The stellar mass function itself has
uncertainty due to effects such as cosmic variance (e.g., Moster et al., 2011;
Steinhardt et al., 2021), SED modeling uncertainties (e.g., Marchesini et al.,
2009; Pacifici et al., 2023), photometric redshift, and Poisson noise. The
weight fluctuates according to the uncertainty of the stellar mass function to
include these effects.

Several studies report that the size-mass relation bends at
$\log{(M_{\star}/M_{\odot})}=10-10.5$ at $z<2$ (e.g., Mowla et al., 2019a;
Nedkova et al., 2021; Kawinwanichakij et al., 2021; Cutler et al., 2022). To
avoid the deviation from our assumption in Equation 2, we only fit the
galaxies with $\log{(M_{\star}/M_{\odot})}>10.3$. This is the same stellar
mass limit imposed in van der Wel et al. (2014). In addition, dusty star-
forming galaxies could be contaminating our sample (e.g., Pérez-González et
al., 2023). They could have larger rest-frame optical sizes than typical
quiescent galaxies (e.g., $\sim 2\,{\rm kpc}$ at $z\sim 3$ in Gillman et al.,
2023), driving a larger scatter of the relation. Thus, we randomly remove 20%
of the galaxies from the sample in the fitting, which is the contamination
fraction of $UVJ$ quiescent galaxies in Schreiber et al. (2018). Following the
literature (e.g., van der Wel et al., 2014), we randomly select those to be
removed by weighting galaxies by the ratio of the stellar mass function of
star-forming galaxies and that of quiescent galaxies at $3<z<3.5$ in Weaver et
al. (2022). Even if we do not account for these potential contaminants, we
find that the best-fit parameters are consistent within their $1\sigma$
uncertainty. In the fitting, the Markov chain Monte Carlo (MCMC) method is
employed using emcee (Foreman-Mackey et al., 2013). The fitting procedure,
including fluctuating the weight $W$ and randomly removing the possible
contaminants, is repeated 500 times. The best-fit values of A (the effective
radius at $M_{\star}=5\times 10^{10}M_{\odot}$), $\alpha$ (the slope of the
size-mass relation), and $\sigma_{\log{R_{\rm eff}}}$ (the intrinsic scatter
of the size-mass relation) are obtained from the total of their chain.

Figure 6 shows the best-fit of the size-mass relation, and Table 1 summarizes
the best-fit values of the parameters. The best-fit $\alpha$ confirms that the
significant positive correlation between the stellar mass and the effective
radius exists in both $UVJ$-selected and $NUVUVJ$-selected galaxies, even
considering the $1\sigma$ uncertainty. The intrinsic scatter is $\sim 0.2$ dex
for both the $UVJ$-selected and $NUVUVJ$-selected sample, similar to those
reported at $z<3$ (e.g., van der Wel et al., 2014; Kawinwanichakij et al.,
2021).

We split $UVJ$-selected quiescent galaxies into two subsamples based on their
photometric redshift, i.e., $z<3.3$ and $z>3.3$. The threshold corresponds to
the median of the redshift of the sample. Figure 7 show the size-mass relation
of these subsamples and their best fit. Table 1 summarizes the best-fit values
of the parameters. Though the uncertainty is more considerable than the total
$UVJ$-selected sample, we see a positive correlation between the stellar mass
and the effective radius even in these subsamples. The effective radius at the
fixed stellar mass is smaller in the high-redshift sample than in the low-
redshift sample. This suggests a potential redshift evolution of the sizes of
quiescent galaxies at $z\geq 3$.

Table 1: Summary of the obtained best-fit parameters of the size-mass relation. Sample | NaaThe number of galaxies used in deriving the best-fit parameters of the size-mass relation. |  $\log{(A/{\rm kpc})}$bbTheir uncertainties are taken from the 16th to 84th percentile of the total chain in the MCMC fitting. |  $\alpha$bbTheir uncertainties are taken from the 16th to 84th percentile of the total chain in the MCMC fitting. | $\sigma_{\log{R_{\rm eff}}}\tablenotemark{b}$
---|---|---|---|---
$UVJ$ | 22 | $-0.18_{-0.06}^{+0.06}$ | $0.32_{-0.21}^{+0.20}$ | $0.18_{-0.05}^{+0.07}$
$UVJ$ ($z<3.3$) | 12 | $-0.14_{-0.08}^{+0.08}$ | $0.43_{-0.35}^{+0.29}$ | $0.16_{-0.06}^{+0.09}$
$UVJ$ ($z>3.3$) | 10 | $-0.26_{-0.17}^{+0.15}$ | $0.39_{-0.49}^{+0.43}$ | $0.28_{-0.11}^{+0.13}$
$NUVUVJ$ | 12 | $-0.24_{-0.09}^{+0.08}$ | $0.54_{-0.28}^{+0.24}$ | $0.20_{-0.08}^{+0.13}$
vDW14 ($z=2.75$)ccThe best-fit parameters of early-type galaxies at $2.5<z<3.0$ reported in van der Wel et al. (2014). | - | $-0.06\pm 0.03$ | $0.79\pm 0.07$ | $0.14\pm 0.03$
Figure 7: Size-mass relation of $UVJ$-selected quiescent galaxies split into
two subsamples, $z<3.3$ (left) and $z>3.3$ (right). The meanings of the
markers and the lines are the same as in Figure 6.

###  4.3 Evolution of the size-mass relation of quiescent galaxies from $z\sim
4$

Figure 8 shows the evolution of the effective radius at $M_{\star}=5\times
10^{10}\,M_{\odot}$ (i.e., $A$) and the slope $\alpha$. In this figure, we
compare our values with those of quiescent galaxies at lower redshift ($z<3$)
(van der Wel et al., 2014; Mowla et al., 2019b; Nedkova et al., 2021). We note
that all of the above studies select quiescent galaxies based on the $UVJ$
diagram. We confirm that the effective radius at $M_{\star}=5\times
10^{10}\,M_{\odot}$ of this study is smaller than those at $z<3$, implying
that the size of the quiescent galaxies monotonically increases since $3<z<4$.
Figure 9 also shows the redshift evolution of the size-mass relation. Our
measured relations are $\sim 0.2$ dex below from those at $z\sim 2$ from the
literature and have the smallest amplitude in the studies. Moreover, our
measurement is in good agreement with the extrapolation of the redshift
evolution of the effective radii suggested in van der Wel et al. (2014), which
is $R_{e}=5.6\times(1+z)^{-1.48}\,{\rm kpc}$ or $R_{e}=4.3\times
h(z)^{-1.29}\,{\rm kpc}$ (see Figure 8). The $UVJ$-selected sample subdivided
by redshift is also consistent with these extrapolations.

The slope of the size-mass relation of quiescent galaxies does not seem to
evolve with redshift. Those drawn from $NUVUVJ$\- and $UVJ$-selected quiescent
galaxies - also divided in redshift bins - are consistent. Our values are also
consistent with those at $z<3$ reported in Mowla et al. (2019b) and Nedkova et
al. (2021). Also, this slope is consistent with the high-mass end slope of the
relation in Cutler et al. (2022). However, our values are slightly flatter
than those in van der Wel et al. (2014), in particular for the $UVJ$-selected
sample. The possible origin of this difference will be discussed in Section 6.

Lastly, Figure 9 compares our size-mass relation with the size measurements of
individual quiescent galaxies at $z\geq 2.5$ in the literature (Gobat et al.,
2012; Kubo et al., 2018; Saracco et al., 2020; Lustig et al., 2021; Esdaile et
al., 2021; Forrest et al., 2022; Carnall et al., 2023b). Using the
$K_{s}$-band image of the AO-assisted VLT/HAWK-I, we also obtain the
structural parameter of a quiescent galaxy at $z=4.01$ (SXDS-27434) reported
in Tanaka et al. (2019) and Valentino et al. (2020). Its effective radius is
estimated as $r_{\rm e}=0.74^{+0.20}_{-0.09}\,{\rm kpc}$ (see Appendix C for
further details). Though the observed wavelength differs among the sample, and
most of these values are not corrected for the wavelength dependence of the
size, they are overall consistent with our measurements. It should be noted
that all $z\geq 4$ measurements in Kubo et al. (2018) and Carnall et al.
(2023b) are below our size-mass relations, so as the SXDS-27434 at $z=4.01$.
This is in agreement with our $z\geq 4$ quiescent galaxies (#2876 at
$z=4.146$, #185 at $z=4.508$, and #10084 at $z=4.633$) with two-thirds smaller
than our relations, which are $0.42\,{\rm kpc}$ at
$\log{(M_{\star}/M_{\odot})}=11.0$ on average. Such a systematically small
size for the $z\geq 4$ sample may imply that the size evolution of quiescent
galaxies starts at $z\geq 4$.

Figure 8: Left panel: Redshift evolution of the effective radius at
$M_{\star}=5\times 10^{10}\,M_{\odot}$ in the best-fit size-mass relation.
Large red and blue circles show the cases for $UVJ$-selected and
$NUVUVJ$-selected quiescent galaxies. Small red circles are the case for
redshift-sub binned $UVJ$-selected quiescent galaxies. Black circles, gray
squares, and open triangles are the lower redshift measurements in van der Wel
et al. (2014), Mowla et al. (2019b), and Nedkova et al. (2021), respectively.
Solid and dotted lines are the size evolution function suggested in van der
Wel et al. (2014), and the dashed lines are their extrapolation at $z>3$.
Right panel: Redshift evolution of the slope of the size-mass relation
($\alpha$). The symbols’ meanings and the data points’ horizontal axis values
are the same as in the left panel. Figure 9: Compilation of the size-mass
relation of quiescent galaxies up to $z\sim 4$. Those from this study ($z\geq
3$) are shown in the solid and dashed lines for $UVJ$ and $NUVUVJ$-selected
samples with the uncertainty with the shaded region, respectively. Those at
lower redshift are taken from the various literature, van der Wel et al.
(2014) in dotted lines ($0<z<3$), Kawinwanichakij et al. (2021) in dashed-
dotted lines ($0<z<1$), and Nedkova et al. (2021) in double-dotted line
($0<z<2$), respectively. Markers represent the size measurement of individual
quiescent galaxies from the literature at $z\geq 2.5$ (Gobat et al., 2012;
Kubo et al., 2018; Saracco et al., 2020; Lustig et al., 2021; Esdaile et al.,
2021; Forrest et al., 2022; Carnall et al., 2023b). We note that none of these
data points are corrected for the wavelength dependence of the size and the
wavelength of the used filter is different among them. If the rest-frame
wavelength of the images used for deriving the size is less than $0.4\,{\rm\mu
m}$, the color of their marker edge and error bar is set to be gray. If it is
longer than $0.4\,{\rm\mu m}$, it is set to be black.

##  5 Morphological Properties

###  5.1 Sérsic Index and Disk Fraction

The Sérsic index measured in F277W, the filter close to rest-frame
$0.5\,{\rm\mu m}$ at $z\sim 3-4$, is shown in the top panel of Figure 10.
Here, we only focus on the objects with successful measurements of the Sérsic
index (i.e., galaxies whose best fits are not located at the parameter range’s
edge). The measured Sérsic indices range from $n=0.5$ to $n=6$, with median
values of $n=3.42_{-0.57}^{+0.49}$ and $n=3.61_{-0.47}^{+0.51}$ for
$UVJ$-selected and $NUVUVJ$-selected quiescent galaxies, respectively. These
values and uncertainties are computed by bootstrapping 5000 times, considering
the uncertainty of the Sérsic index of each galaxy. The top panel of Figure 10
shows the median of two subsamples divided according to their redshift, as
done in Section 4.2. They do not show significant redshift evolution in their
median values and are in good agreement with lower redshift results (e.g., van
der Wel et al., 2012; Marsan et al., 2019; Lustig et al., 2021; Esdaile et
al., 2021). This suggests that most quiescent galaxies have bulge-dominant
structures as early as $z\sim 4$.

One point that we should mention from our sample is that some quiescent
galaxies have a low Sérsic index as $n\sim 0.5-2$. To investigate this
further, we next derive the disk fraction, defined as the fraction of objects
with $n<2$, i.e., the number of the galaxies with $n<2$ divided by the total
number of the galaxies. These are also computed by bootstrapping. The bottom
panel of Figure 10 shows the disk fraction as a function of redshift. Though
both bins at lower redshift are consistent with zero, the disk fraction
increases towards higher redshift. In particular, the high redshift bin
($z>3.3$) for $UVJ$-selected quiescent galaxies has a non-zero value
($0.36_{-0.18}^{+0.09}$) even considering the $1\sigma$ uncertainty. Our
sample suggests that some fraction of quiescent galaxies at high redshift can
be disk-dominant morphology, in line with the literature (e.g., Toft et al.,
2017; Newman et al., 2018; Fudamoto et al., 2022).

Figure 10: Top panel: Sérsic index in F277W filter as a function of the
redshift. The small colored circles are the individual quiescent galaxies
selected by $UVJ$ or $NUVUVJ$. The larger circles and squares are the median
of them at $z<3.3$ and $z>3.3$. Black markers show the median of the
measurement of quiescent galaxies at $z<3$ in the literature (van der Wel et
al., 2012; Marsan et al., 2019; Lustig et al., 2021; Esdaile et al., 2021).
The stellar mass range of their sample is $\log{(M_{\star}/M_{\odot})}>10.3$,
$\log{(M_{\star}/M_{\odot})}>11.25$, $10.8<\log{(M_{\star}/M_{\odot})}<11.3$,
and $11.0<\log{(M_{\star}/M_{\odot})}<11.31$ for van der Wel et al. (2012),
Marsan et al. (2019), Lustig et al. (2021), and Esdaile et al. (2021),
respectively. Quiescent galaxies in van der Wel et al. (2012) are selected in
the same manner as van der Wel et al. (2014). Their values are computed by
bootstrapping, the same as the value from this work. Bottom panel: Disk
fraction, defined as the fraction of sources with F277W Sérsic index below 2.
Circles and squares represent $UVJ$ and $NUVUVJ$-selected quiescent galaxies,
respectively.

###  5.2 Axis Ratio

The observed axis ratio of each galaxy projected to the images depends on its
inclination angle, but by taking the median and removing the effect from the
inclination angle, we trace the intrinsic shape of the sample. The median axis
ratio in F277W is $q=0.61_{-0.09}^{+0.12}$ and $q=0.78_{-0.20}^{+0.04}$ for
$UVJ$-selected and $NUVUVJ$-selected sample, respectively. They are again
obtained from bootstrapping, as for the median value of the Sérsic index
(Section 5.1). The $NUVUVJ$-selected sample has a higher median value than
that of the $UVJ$-selected sample, though with a notably larger uncertainty.
Figure 11 shows the individual values of our sample as a function of their
redshift and the median value for $2.8<z\leq 3.3$ and $3.3<z<4.6$ sample. We
do not see any significant dependence of the axis ratio on either the redshift
or the stellar mass.

Our median values are in good agreement with the measurement of quiescent
galaxies at $z<3$ (e.g., van der Wel et al., 2012; Straatman et al., 2015;
Lustig et al., 2021). At $z\geq 3$, Esdaile et al. (2021) derive the axis
ratio of four quiescent galaxies based on the HST/F160W images. Their median
($q\sim 0.5-0.7$) is also consistent with ours. Coupled with the median Sérsic
index, these values provide further evidence that the majority of quiescent
galaxies have bulge-like structures as early as $z\sim 4$.

Figure 11: Axis ratio in F277W filter as a function of the redshift. Small
colored circles are those of the quiescent galaxies used in this study, and
the larger circles and squares are the median of the $UVJ$-selected and
$NUVUVJ$-selected samples, respectively. They are subdivided into two samples
according to their redshift ($z\leq 3.3$ and $z>3.3$). Black markers show the
measurement of quiescent galaxies in the literature (van der Wel et al., 2012;
Straatman et al., 2015; Marsan et al., 2019; Lustig et al., 2021; Esdaile et
al., 2021).

##  6 Discussion and Summary

This study investigates the morphological properties of 26 quiescent galaxies
at $3<z<5$ in public $JWST$ fields. Quiescent galaxies are identified by
Valentino et al. (2023), selecting targets from all available public $JWST$
data in Cycle 1 based on two rest-frame color selections, i.e., the
traditional $UVJ$ diagram (Williams et al., 2009) and a novel selection from
the rest-frame $NUV$, $U$, $V$, and $J$ band magnitude following Gould et al.
(2023). By applying the Galfit software to all NIRCam images of these
quiescent galaxies, we derive their effective radius, the Sérsic index, and
the axis ratio at observed $0.9-4.4\,{\rm\mu m}$.

Thanks to the measurements in the multiple bands, we derive the rest-frame
$0.5\,{\rm\mu m}$ effective radius by considering the wavelength dependence of
the effective radius. As those at lower redshift (e.g., van der Wel et al.,
2014; Kawinwanichakij et al., 2021; Suess et al., 2022), quiescent galaxies at
$z\geq 3$ show a significant negative color gradient in effective radius,
i.e., having a more compact size at a longer wavelength. This can be due to
the radial dependence of their stellar population. If the inside-out quenching
happens, the older stellar population is more concentrated than younger
stellar populations. Since longer wavelength (i.e., $\lambda>0.5\ \mu$m in the
rest-frame in our case) sees a more evolved stellar population, this can cause
the observed wavelength dependence of the size. Our observed trend implies
that correcting for the wavelength dependence of the size measurements is
critical for fairly comparing the sizes of quiescent galaxies at different
redshifts.

We derive the relation of the stellar mass and the effective radius at the
rest-frame $0.5\,{\rm\mu m}$ for quiescent galaxies at $3<z<5$ and the best-
fit parameters assuming a power-law relation in the stellar mass range of
$\log{(M_{\star}/M_{\odot})}>10.3$. These measurements suggest a significant
correlation between the stellar mass and effective radii as early as $z\geq 3$
with the power law slope of $0.3-0.5$. In addition, by splitting the
$UVJ$-selected sample into two redshift subsamples, we see the redshift
evolution of its amplitude at $z\geq 3$.

The effective radii at $M_{\star}=5\times 10^{10}\,M_{\odot}$ from the best-
fit are $\sim 0.6\,{\rm kpc}$, $\sim 0.2$ dex smaller than that of $z\sim 2$.
This value is the smallest at $0<z<4$. Therefore, we show that the size
evolution of quiescent galaxies has started already at $3<z<5$. Our results
agree with the redshift evolution in size extrapolated from lower redshift
measurements in van der Wel et al. (2014).

The slope of the size-mass relation does not show the significant redshift
evolution in our sample. Though having considerable uncertainty, our
measurement shows a flatter slope ($\sim 0.3-0.5$) than that of van der Wel et
al. (2014, $\sim 0.7-0.8$) in particular for the $UVJ$-selected sample, and
close to those of Mowla et al. (2019b), Nedkova et al. (2021) and the high-
mass end slope of the broken power law fit in Cutler et al. (2022). There can
be several explanations for the possible difference between our results and
van der Wel et al. (2014). Firstly, the method of collecting the wavelength
dependence of the size can make the relation different. We correct it for
individual galaxies so as Nedkova et al. (2021), whereas van der Wel et al.
(2014) apply the same correction factor for all galaxies. If the amount of the
wavelength dependence of the size is different among galaxies, such as if
there is a correlation between the stellar mass and the size gradient like at
lower redshift Kawinwanichakij et al. (2021) or if there is a scatter (see
Section 3.5), this could make the relation in van der Wel et al. (2014)
different from ours. Also, this can be due to the increasing fraction of
transitional galaxies from star-forming to quiescent, in particular at lower
stellar mass. For example, Suess et al. (2021) shows that the median size of
green valley galaxies lies between star-forming and quiescent galaxies. If
this trend persists at $z\geq 3$, such larger galaxies flatten the size-mass
relation. Lastly, this can be due to the possible redshift dependence of the
stellar mass distribution of the sample, i.e., less low-mass galaxies at
higher redshift (see Figure 7). Since the stellar mass limit is well below the
focused stellar mass range (see the discussion in Valentino et al., 2023),
this is likely due to the small number of the sample. Even though this can not
fully explain the difference between our results and van der Wel et al. (2014)
since the slope is still flatter when the sample is split into two redshift
subsamples (Figure 7), this may increase the size at the low stellar mass
compared to high mass due to the larger galaxies at lower redshift, leading to
the flatter slope.

The median Sérsic index and the axis ratios are $n=3.42_{-0.57}^{+0.49}$
($n=3.61_{-0.47}^{+0.51}$) and $q=0.61_{-0.09}^{+0.12}$
($q=0.78_{-0.20}^{+0.04}$) for $UVJ$ ($NUVUVJ$)-selected sample, respectively.
This high Sérsic index and high axis ratio imply that most quiescent galaxies
have bulge-like morphologies. At the same time, we see an increasing trend of
the fraction of disk-like morphology, especially in the $UVJ$-selected sample.
The fact that one of the two spectroscopically confirmed quiescent galaxies at
$z>4$ so far (Tanaka et al., 2019; Valentino et al., 2020; Carnall et al.,
2023b), SXDS-27434, has a low Sérsic index ($n={1.23}^{+0.65}_{-0.91}$, see
Appendix C) also agrees with this trend. These might imply that $3<z<5$ is the
epoch when some fraction of quiescent galaxies is in their morphological
transformation phase from disk galaxies to elliptical galaxies.

This study demonstrates that $JWST$ enables us to explore the structural
properties of quiescent galaxies $z>3$. This study serves as a first look at
this topic with $JWST$. We only focus on the half-light radii here, but note
that future work based on the $JWST$ multi-band imaging will enable us to
determine the half-mass radii, as done with $HST$ (Suess et al., 2021). In
addition, the wider coverage incoming from new imaging surveys with JWST, such
as COSMOS-Web (Casey et al., 2022), will constrain the structural properties
of a significantly larger sample of quiescent galaxies at $z\geq 3$. This will
reduce the uncertainty of the analytical fitting of the size-mass relation.
Also, the systematic comparison with the morphological properties of star-
forming galaxies at $z>3$ will be helpful to constrain their quenching
mechanism, as done at $z<3$ (e.g., Faisst et al., 2017). Moreover, by
increasing the number of quiescent galaxies at $z\geq 4$, we can constrain the
trend at an even higher redshift. In particular, exploring the evolution of
the size-mass relation and the disk fraction will provide a statistical census
into the connection between structural transformation and quenching of the
current most distant quiescent galaxies.

This study is based on the observations associated with programs ERS #1324,
1345, and 1355; ERO #2736; GO #1837 and 2822; GTO #2738; and COM #1063. The
authors acknowledge the teams and PIs for developing their observing program
with a zero-exclusive-access period. Part of the study is based on
observations collected at the European Southern Observatory under ESO program
0104.B-0213(A). This study was supported by JSPS KAKENHI Grant Numbers
JP21K03622, JP22J00495, and JP23K13141. OI acknowledges the funding of the
French Agence Nationale de la Recherche for the project iMAGE (grant
ANR-22-CE31-0007). The Cosmic Dawn Center of Excellence is funded by the
Danish National Research Foundation under grant No. 140. GEM acknowledges the
Villum Fonden research grant 13160 “Gas to stars, stars to dust: tracing star
formation across cosmic time,” grant 37440, “The Hidden Cosmos”.

##  Appendix A Selection of stars for smoothing PSFs and an example of
smoothed PSF

Stars used to compare with the model are selected based on the following four
criteria. (1). Objects with the signal-to-noise ratio of $S/N>200$ in
$0.5\arcsec$ aperture magnitude are selected to only focus on those with
reliable morphology information. We note that saturated stars are not selected
here. (2). Ellipticity, defined as the $(a-b)/(a+b)$, where $a,\ {\rm and}\ b$
are the major and minor axis radius, should be lower than 0.1 to remove the
galaxies. (3). Objects on a sequence of unresolved objects are selected based
on FLUX_RADIUS from sep. We define the range of the sequence as the $1\sigma$
from the mean of the sequence. (4). If objects have brighter sources within
1″(2″) in short (long) wavelength images, they are not used to avoid the
impact of the surrounding sources.

Figure 12 shows the smoothed PSFs used in this study, the natural stars, and
the original PSFs from WebbPSF of CEERS as an example. The original PSF from
WebbPSF is clearly more compact than the natural stars, particularly at
shorter wavelengths, implying the need for smoothing. The PSFs smoothed by the
Moffat function agree mostly well with natural stars in all filters.

Figure 12: Summary of point spread function of F115W, F150W, F200W, F277W,
F356W, F410M and F444W filter in CEERS. Each panel has the radial profile of
our smoothed PSF from WebbPSF (red), the original PSF from WebbPSF (blue), the
median-stacked stars (black), and individual stars (thin gray) at the top. All
of them are normalized to have the unity flux at $r=0\arcsec$. The bottom
panel shows the ratio of the radial profile between the smoothed PSF (original
PSF) and the stacked stars in the red (blue) line. The ratio between those and
the individual stars is shown in thinner lines.

##  Appendix B Summary of fitting result

Figure 13, Figure 14, and Figure 15 show the fitting results in F200W and
F277W of all 26 objects which are used in this study. In addition, Table 2
summarizes the best-fit parameters obtained in fitting F277W images and the
effective radius at rest-frame $0.5\,{\rm\mu m}$ and the wavelength dependence
of the size derived in Section 3.5.

Figure 13: Summary of the fitting result in F200W (left) and F277W (right)
images, respectively. For each filter and each object, its image, mask, model,
and residual are shown. All are $4.0\arcsec\times 4.0\arcsec$. The images are
scaled by their noise ($\pm 15\sigma$). The black contours start at $10\%$ of
the peak pixel count of each object and increase by $\sqrt{2}$. Figure 14:
Continued from Figure 13. Figure 15: Continued from Figure 13. Table 2:
Summary of Galfit fitting.

IDaaThe IDs are identical to those in Valentino et al. (2023). | field |  $z_{\rm phot}$ bbThey are derived with eazy-py in Valentino et al. (2023). |  $\log{(M_{\star}/M_{\odot})}$bbThey are derived with eazy-py in Valentino et al. (2023). | $r_{\rm eff,f277w}$ | $n_{\rm f277w}$ | $q_{\rm f277w}$ | $PA_{\rm f277w}$ | $r_{\rm eff,0.5\mu m}\tablenotemark{c}$ |  $\gamma$ddThe wavelength dependence of the effective radius derived by Equation 1. | flageefootnotemark:
---|---|---|---|---|---|---|---|---|---|---
|  |  |  | (kpc) |  |  | (deg) | (kpc) |  |
7174 | CEERS | 2.83 | 10.34 | $0.621\pm 0.021$ | $2.25\pm 0.23$ | $0.701\pm 0.025$ | $-61\pm 3$ | $0.798\pm 0.057$ | $-0.57\pm 0.11$ | 0
19471 | J1235 | 2.83 | 10.46 | $0.269\pm 0.006$ | $5.65\pm 0.4$ | $0.315\pm 0.012$ | $14\pm 1$ | $0.464\pm 0.064$ | $-0.87\pm 0.3$ | 0
788 | SGAS1723 | 2.85 | 10.72 | $3.744\pm 0.988$ | $4.57\pm 0.96$ | $0.365\pm 0.034$ | $31\pm 3$ | $4.407\pm 3.339$ | $-0.73\pm 1.02$ | 0
2278 | PRIMER | 2.85 | 10.41 | $0.722\pm 0.094$ | $5.52\pm 0.97$ | $0.844\pm 0.055$ | $-78\pm 15$ | $0.893\pm 0.279$ | $-0.49\pm 0.44$ | 0
4675 | CEERS | 2.91 | 10.78 | $0.5\pm 0.012$ | $4.78\pm 0.35$ | $0.514\pm 0.01$ | $73\pm 1$ | $0.602\pm 0.032$ | $-0.51\pm 0.08$ | 0
2809 | PRIMER | 2.92 | 10.77 | $0.823\pm 0.035$ | $2.89\pm 0.27$ | $0.508\pm 0.016$ | $77\pm 1$ | $0.926\pm 0.066$ | $-0.3\pm 0.1$ | 0
28463 | QUINTET | 2.93 | 10.4 | $1.023\pm 0.2$ | $[4.0]$ | $0.626\pm 0.08$ | $-32\pm 8$ | $0.806\pm 0.402$ | $0.82\pm 0.74$ | 2
2560 | PRIMER | 2.93 | 11.24 | $1.587\pm 0.037$ | $2.25\pm 0.08$ | $0.843\pm 0.01$ | $-45\pm 2$ | $1.836\pm 0.113$ | $-0.34\pm 0.09$ | 0
7399 | CEERS | 2.98 | 10.21 | $2.508\pm 0.202$ | $[4.0]$ | $0.243\pm 0.02$ | $37\pm 1$ | $4.882\pm 1.185$ | $-1.74\pm 0.38$ | 2
277 | PRIMER | 3.12 | 10.62 | $0.336\pm 0.012$ | $2.11\pm 0.31$ | $0.86\pm 0.043$ | $35\pm 11$ | $0.35\pm 0.021$ | $-0.07\pm 0.11$ | 0
2796 | NEP-TDF | 3.13 | 10.53 | $0.809\pm 0.047$ | $3.88\pm 0.35$ | $0.817\pm 0.02$ | $75\pm 3$ | $0.952\pm 0.103$ | $-0.65\pm 0.17$ | 0
11266 | CEERS | 3.18 | 10.74 | $0.494\pm 0.018$ | $4.89\pm 0.41$ | $0.862\pm 0.017$ | $61\pm 4$ | $0.575\pm 0.046$ | $-0.44\pm 0.13$ | 0
7955 | CEERS | 3.21 | 10.54 | $0.517\pm 0.011$ | $3.52\pm 0.26$ | $0.566\pm 0.012$ | $77\pm 1$ | $0.591\pm 0.046$ | $-0.57\pm 0.14$ | 0
9621 | CEERS | 3.38 | 10.98 | $1.272\pm 0.015$ | $1.7\pm 0.05$ | $0.421\pm 0.006$ | $9\pm 0$ | $1.374\pm 0.028$ | $-0.37\pm 0.04$ | 0
15782 | CEERS | 3.35 | 11.39 | $1.711\pm 0.043$ | $4.05\pm 0.1$ | $0.806\pm 0.005$ | $79\pm 1$ | $1.857\pm 0.073$ | $-0.38\pm 0.07$ | 0
35943 | QUINTET | 3.36 | 10.12 | $0.708\pm 0.067$ | $1.98\pm 0.43$ | $0.905\pm 0.063$ | $-81\pm 82$ | $0.719\pm 0.134$ | $-0.11\pm 0.32$ | 0
8083 | CEERS | 3.41 | 10.87 | $0.672\pm 0.013$ | $3.17\pm 0.13$ | $0.775\pm 0.009$ | $-71\pm 1$ | $0.719\pm 0.023$ | $-0.31\pm 0.06$ | 0
1863 | PRIMER | 3.46 | 10.72 | $0.739\pm 0.043$ | $2.58\pm 0.38$ | $0.447\pm 0.029$ | $-28\pm 2$ | $0.848\pm 0.126$ | $-0.41\pm 0.26$ | 1
8967 | NEP-TDF | 3.47 | 9.77 | $2.214\pm 0.153$ | $0.52\pm 0.1$ | $0.502\pm 0.035$ | $60\pm 3$ | $2.159\pm 0.635$ | $-0.37\pm 0.56$ | 0
8782 | CEERS | 3.53 | 10.39 | $0.316\pm 0.012$ | $4.85\pm 0.53$ | $0.941\pm 0.028$ | $46\pm 16$ | $0.332\pm 0.031$ | $-0.18\pm 0.18$ | 0
9622 | CEERS | 3.54 | 11.11 | $0.487\pm 0.005$ | $3.64\pm 0.14$ | $0.485\pm 0.006$ | $60\pm 0$ | $0.532\pm 0.031$ | $-0.11\pm 0.12$ | 0
1427 | PRIMER | 3.55 | 10.18 | $0.086\pm 0.016$ | $[4.0]$ | $0.514\pm 0.205$ | $-37\pm 13$ | $0.054\pm 0.024$ | $0.67\pm 0.95$ | 3
10378 | CEERS | 3.67 | 10.83 | $1.002\pm 0.081$ | $5.34\pm 0.53$ | $0.648\pm 0.016$ | $70\pm 2$ | $1.161\pm 0.106$ | $-0.85\pm 0.17$ | 0
2876 | CEERS | 4.15 | 10.81 | $0.722\pm 0.014$ | $1.67\pm 0.11$ | $0.464\pm 0.016$ | $17\pm 1$ | $0.743\pm 0.019$ | $-0.45\pm 0.06$ | 0
185 | PRIMER | 4.51 | 10.91 | $0.033\pm 0.035$ | $4.66\pm 3.52$ | $0.232\pm 0.313$ | $-82\pm 81$ | $0.034\pm 0.046$ | $0\pm 4\times 10^{7}$ | 1
10084 | CEERS | 4.63 | 11.23 | $0.579\pm 0.016$ | $[4.0]$ | $0.968\pm 0.015$ | $34\pm 16$ | $0.475\pm 0.038$ | $0.15\pm 0.21$ | 2

eefootnotemark:

The quality of fitting. flag$=0$: good fit. flag$=1$: The Galfit reports an
uncertain fit for at least one parameter in one band (i.e., it has parameters
marked by *…* in the output file of Galfit). flag$=2$: Its Sérsic index is
fixed as $n=4.0$. flag$=3$: It is flagged as both 1 and 2.

##  Appendix C New size estimate of SXDS-27434 at $z=4.01$

Here we report the size estimate at rest-frame $\sim 4200$ Å of the galaxy
SXDS-27434 at $z=4.01$ (Tanaka et al., 2019; Valentino et al., 2020). This is
based on observations with the High Acuity Wide-field K-band Imager (HAWK-I,
Pirard et al., 2004; Kissler-Patig et al., 2008) coupled with the Ground Layer
Adaptive Optic system GRAAL (Paufique et al., 2010) (Program ID 0104.B-0213,
PI: F. Valentino). The target ($\mathrm{R.A.,Dec.}=34.29871,-4.98987$) was
observed with the $K_{\rm s}$ filter for $2040\,\mathrm{s}$ on source in
service mode on October 22nd, 2019. We applied the standard jitter scheme
within a 25″ box and obtained 34 offset exposures ($\mathrm{NDIT\times
DIT}=6\times 10\,\mathrm{s}$). We reduce the data as described in Brammer et
al. (2016)222https://github.com/gbrammer/HAWKI-FF. Notably, particular care is
given to an optimal sky subtraction with a “sky flat” obtained from the median
of all the science frames. Finally, the images are “drizzled” (Fruchter &
Hook, 2002) to a pixel size of $0\farcs 08$ and their flux calibration tied to
the catalog described in Tanaka et al. (2019). We model the surface brightness
of the target as a single Sérsic profile with galfit. We create a natural PSF
by stacking unsaturated stars falling in the HAWK-I footprint, selected in the
stellar locus in the FLUX_RADIUS-MAG_AUTO plane generated with Source
Extractor (Bertin & Arnouts, 1996) and visually inspected. The final average
seeing is of $\mathrm{FWHM}=0\farcs 34$. The best-fit model and the residuals
are shown in Figure 16. We measure an effective semi-major axis of $r_{\rm
e}=0.74^{+0.20}_{-0.09}$ kpc and a Sérsic index $n={1.23}^{+0.65}_{-0.91}$
with an axis ratio $q=0.88^{+0.06}_{-0.26}$. The uncertainties are estimated
by bootstrapping $1000\times$ the input cutouts within the rms per pixel
measured from the calibrated image and relaunching galfit in the same
configuration as for the science measurement. The error bars are computed as
the 16-84% percentiles of the distribution of the parameters in the
bootstrapped images. The measurements are robust against the choice of the
size of the drizzled pixels ($0\farcs 0.8$ or $0\farcs 1$) and a fixed or
freely varying sky level. The target is resolved, and the radius well-
constrained. This is confirmed by the structured residuals obtained when
fitting a single point source in lieu of a Sérsic profile. The effective
radius estimate is consistent with that from $i$-band Hyper-SuprimeCam
observations reported in Tanaka et al. (2019) ($r_{\rm e}=0.76\pm 0.20$ kpc).

Figure 16: HAWK-I+GRAAL image (left), galfit best-fit Sérsic model (center),
and residuals (right) of SXDS-27434 at $z=4.01$. The size of each cutout is
10″. The average seeing size is shown by the white circle in the left panel.

## References

  * Bagley et al. (2023) Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJ, 946, L12, doi: 10.3847/2041-8213/acbb08
  * Barbary (2016) Barbary, K. 2016, The Journal of Open Source Software, 1, 58, doi: 10.21105/joss.00058
  * Barro et al. (2013) Barro, G., Faber, S. M., Pérez-González, P. G., et al. 2013, ApJ, 765, 104, doi: 10.1088/0004-637X/765/2/104
  * Bertin (2011) Bertin, E. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 442, Astronomical Data Analysis Software and Systems XX, ed. I. N. Evans, A. Accomazzi, D. J. Mink, & A. H. Rots, 435
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * 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
  * Bradley et al. (2022) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, astropy/photutils: 1.5.0, 1.5.0, Zenodo, doi: 10.5281/zenodo.6825092
  * Brammer (2021) Brammer, G. 2021, eazy-py, 0.5.2, doi: 10.5281/zenodo.5012704
  * Brammer & Matharu (2021) Brammer, G., & Matharu, J. 2021, gbrammer/grizli: Release 2021, 1.3.2, Zenodo, Zenodo, doi: 10.5281/zenodo.1146904
  * Brammer et al. (2022) Brammer, G., Strait, V., Matharu, J., & Momcheva, I. 2022, grizli, 1.5.0, Zenodo, Zenodo, doi: 10.5281/zenodo.6672538
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786
  * Brammer et al. (2016) Brammer, G. B., Marchesini, D., Labbé, I., et al. 2016, ApJS, 226, 6, doi: 10.3847/0067-0049/226/1/6
  * 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
  * Casey et al. (2022) Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2022, arXiv e-prints, arXiv:2211.07865, doi: 10.48550/arXiv.2211.07865
  * 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
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * 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
  * D’Eugenio et al. (2021) D’Eugenio, C., Daddi, E., Gobat, R., et al. 2021, A&A, 653, A32, doi: 10.1051/0004-6361/202040067
  * Ding et al. (2022) Ding, X., Silverman, J. D., & Onoue, M. 2022, ApJ, 939, L28, doi: 10.3847/2041-8213/ac9c02
  * Esdaile et al. (2021) Esdaile, J., Glazebrook, K., Labbé, I., et al. 2021, ApJ, 908, L35, doi: 10.3847/2041-8213/abe11e
  * Faisst et al. (2017) Faisst, A. L., Carollo, C. M., Capak, P. L., et al. 2017, ApJ, 839, 71, doi: 10.3847/1538-4357/aa697a
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Forrest et al. (2020) Forrest, B., Marsan, Z. C., Annunziatella, M., et al. 2020, ApJ, 903, 47, doi: 10.3847/1538-4357/abb819
  * Forrest et al. (2022) Forrest, B., Wilson, G., Muzzin, A., et al. 2022, ApJ, 938, 109, doi: 10.3847/1538-4357/ac8747
  * Fruchter & Hook (2002) Fruchter, A. S., & Hook, R. N. 2002, PASP, 114, 144, doi: 10.1086/338393
  * Fudamoto et al. (2022) Fudamoto, Y., Inoue, A. K., & Sugahara, Y. 2022, ApJ, 938, L24, doi: 10.3847/2041-8213/ac982b
  * Gillman et al. (2023) Gillman, S., Gullberg, B., Brammer, G., et al. 2023, arXiv e-prints, arXiv:2303.17246, doi: 10.48550/arXiv.2303.17246
  * Glazebrook et al. (2017) Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Nature, 544, 71, doi: 10.1038/nature21680
  * Gobat et al. (2012) Gobat, R., Strazzullo, V., Daddi, E., et al. 2012, ApJ, 759, L44, doi: 10.1088/2041-8205/759/2/L44
  * Gómez-Guijarro et al. (2023) Gómez-Guijarro, C., Magnelli, B., Elbaz, D., et al. 2023, arXiv e-prints, arXiv:2304.08517, doi: 10.48550/arXiv.2304.08517
  * Gould (2023) Gould, K. M. L. 2023, GMM-quiescent, 1.0.0, doi: 10.5281/zenodo.7712937
  * Gould et al. (2023) Gould, K. M. L., Brammer, G., Valentino, F., et al. 2023, arXiv e-prints, arXiv:2302.10934, doi: 10.48550/arXiv.2302.10934
  * 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
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Kawinwanichakij et al. (2021) Kawinwanichakij, L., Silverman, J. D., Ding, X., et al. 2021, ApJ, 921, 38, doi: 10.3847/1538-4357/ac1f21
  * Kissler-Patig et al. (2008) Kissler-Patig, M., Pirard, J. F., Casali, M., et al. 2008, A&A, 491, 941, doi: 10.1051/0004-6361:200809910
  * Kocevski et al. (2018) Kocevski, D. D., Hasinger, G., Brightman, M., et al. 2018, ApJS, 236, 48, doi: 10.3847/1538-4365/aab9b4
  * Kubo et al. (2018) Kubo, M., Tanaka, M., Yabe, K., et al. 2018, ApJ, 867, 1, doi: 10.3847/1538-4357/aae3e8
  * Lam et al. (2015) Lam, S. K., Pitrou, A., & Seibert, S. 2015, in Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC - LLVM ’15 (New York, New York, USA: ACM Press), 1–6, doi: 10.1145/2833157.2833162
  * Lustig et al. (2021) Lustig, P., Strazzullo, V., D’Eugenio, C., et al. 2021, MNRAS, 501, 2659, doi: 10.1093/mnras/staa3766
  * Marchesini et al. (2009) Marchesini, D., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2009, ApJ, 701, 1765, doi: 10.1088/0004-637X/701/2/1765
  * Marsan et al. (2019) Marsan, Z. C., Marchesini, D., Muzzin, A., et al. 2019, ApJ, 871, 201, doi: 10.3847/1538-4357/aaf808
  * McKinney (2010) McKinney, W. 2010, Data Structures for Statistical Computing in Python, Tech. rep., doi: 10.25080/majora-92bf1922-00a
  * Mosleh et al. (2017) Mosleh, M., Tacchella, S., Renzini, A., et al. 2017, ApJ, 837, 2, doi: 10.3847/1538-4357/aa5f14
  * Moster et al. (2011) Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113, doi: 10.1088/0004-637X/731/2/113
  * Mowla et al. (2019a) Mowla, L., van der Wel, A., van Dokkum, P., & Miller, T. B. 2019a, ApJ, 872, L13, doi: 10.3847/2041-8213/ab0379
  * Mowla et al. (2019b) Mowla, L. A., van Dokkum, P., Brammer, G. B., et al. 2019b, ApJ, 880, 57, doi: 10.3847/1538-4357/ab290a
  * 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
  * Nandra et al. (2015) Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, ApJS, 220, 10, doi: 10.1088/0067-0049/220/1/10
  * Nedkova et al. (2021) Nedkova, K. V., Häußler, B., Marchesini, D., et al. 2021, MNRAS, 506, 928, doi: 10.1093/mnras/stab1744
  * Newman et al. (2018) Newman, A. B., Belli, S., Ellis, R. S., & Patel, S. G. 2018, ApJ, 862, 126, doi: 10.3847/1538-4357/aacd4f
  * Newville et al. (2014) Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, 0.8.0, Zenodo, doi: 10.5281/zenodo.11813
  * Nipoti et al. (2012) Nipoti, C., Treu, T., Leauthaud, A., et al. 2012, MNRAS, 422, 1714, doi: 10.1111/j.1365-2966.2012.20749.x
  * Oke & Gunn (1983) Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713, doi: 10.1086/160817
  * Ono et al. (2022) Ono, Y., Harikane, Y., Ouchi, M., et al. 2022, arXiv e-prints, arXiv:2208.13582, doi: 10.48550/arXiv.2208.13582
  * 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
  * Paufique et al. (2010) Paufique, J., Bruton, A., Glindemann, A., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7736, Adaptive Optics Systems II, ed. B. L. Ellerbroek, M. Hart, N. Hubin, & P. L. Wizinowich, 77361P, doi: 10.1117/12.858261
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266, doi: 10.1086/340952
  * Peng et al. (2010) —. 2010, AJ, 139, 2097, doi: 10.1088/0004-6256/139/6/2097
  * Pérez-González et al. (2023) Pérez-González, P. G., Barro, G., Annunziatella, M., et al. 2023, ApJ, 946, L16, doi: 10.3847/2041-8213/acb3a5
  * Perrin et al. (2014) Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, 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. (2012) Perrin, M. D., Soummer, R., Elliott, E. M., Lallo, M. D., & Sivaramakrishnan, A. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, ed. M. C. Clampin, G. G. Fazio, H. A. MacEwen, & J. Oschmann, Jacobus M., 84423D, doi: 10.1117/12.925230
  * Pirard et al. (2004) Pirard, J.-F., Kissler-Patig, M., Moorwood, A., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5492, Ground-based Instrumentation for Astronomy, ed. A. F. M. Moorwood & M. Iye, 1763–1772, doi: 10.1117/12.578293
  * Pontoppidan et al. (2022) Pontoppidan, K. M., Barrientes, J., Blome, C., et al. 2022, ApJ, 936, L14, doi: 10.3847/2041-8213/ac8a4e
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Robitaille et al. (2013) Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Saracco et al. (2020) Saracco, P., Marchesini, D., La Barbera, F., et al. 2020, ApJ, 905, 40, doi: 10.3847/1538-4357/abc7c4
  * Schreiber et al. (2018) Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018, A&A, 618, A85, doi: 10.1051/0004-6361/201833070
  * Sérsic (1963) Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * 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
  * Steinhardt et al. (2021) Steinhardt, C. L., Jespersen, C. K., & Linzer, N. B. 2021, ApJ, 923, 8, doi: 10.3847/1538-4357/ac2a2f
  * 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. (2019) Suess, K. A., Kriek, M., Price, S. H., & Barro, G. 2019, ApJ, 877, 103, doi: 10.3847/1538-4357/ab1bda
  * Suess et al. (2021) —. 2021, ApJ, 915, 87, doi: 10.3847/1538-4357/abf1e4
  * Suess et al. (2022) Suess, K. A., Bezanson, R., Nelson, E. J., et al. 2022, ApJ, 937, L33, doi: 10.3847/2041-8213/ac8e06
  * Tanaka et al. (2019) Tanaka, M., Valentino, F., Toft, S., et al. 2019, ApJ, 885, L34, doi: 10.3847/2041-8213/ab4ff3
  * Toft et al. (2017) Toft, S., Zabl, J., Richard, J., et al. 2017, Nature, 546, 510, doi: 10.1038/nature22388
  * Tzanavaris et al. (2014) Tzanavaris, P., Gallagher, S. C., Hornschemeier, A. E., et al. 2014, ApJS, 212, 9, doi: 10.1088/0067-0049/212/1/9
  * Valentino (2023) Valentino, F. 2023, ”An atlas of color-selected quiescent galaxies at z¿3 in public JWST fields” (Supplementary material), 1.0, Zenodo, doi: 10.5281/zenodo.7614908
  * Valentino et al. (2020) Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93, doi: 10.3847/1538-4357/ab64dc
  * Valentino et al. (2023) Valentino, F., Brammer, G., Gould, K. M. L., et al. 2023, ApJ, 947, 20, doi: 10.3847/1538-4357/acbefa
  * 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 Dokkum et al. (2015) van Dokkum, P. G., Nelson, E. J., Franx, M., et al. 2015, ApJ, 813, 23, doi: 10.1088/0004-637X/813/1/23
  * 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
  * Weaver et al. (2022) Weaver, J. R., Davidzon, I., Toft, S., et al. 2022, arXiv e-prints, arXiv:2212.02512, doi: 10.48550/arXiv.2212.02512
  * 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
  * Yang et al. (2022) Yang, L., Morishita, T., Leethochawalit, N., et al. 2022, ApJ, 938, L17, doi: 10.3847/2041-8213/ac8803
  * Zhuang & Shen (2023) Zhuang, M.-Y., & Shen, Y. 2023, arXiv e-prints, arXiv:2304.13776, doi: 10.48550/arXiv.2304.13776