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

Kei Ito, Francesco Valentino, Gabriel Brammer, Andreas L. Faisst, Steven Gillman, Carlos Gómez-Guijarro, Katriona M. L. Gould, Kasper E. Heintz, Olivier Ilbert, Christian Kragh Jespersen, et al.

► To cite this version:

Kei Ito, Francesco Valentino, Gabriel Brammer, Andreas L. Faisst, Steven Gillman, et al. Size-Stellar Mass Relation and Morphology of Quiescent Galaxies at $z \ge 3$ in Public JWST Fields. The Astrophysical Journal, 2024, 964, 10.3847/1538-4357/ad2512. insu-04726526

HAL Id: insu-04726526

https://insu.hal.science/insu-04726526v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Kei Ito<sup>1,17</sup> ID, Francesco Valentino<sup>2,3</sup> ID, Gabriel Brammer<sup>3,4</sup> ID, Andreas L. Faisst<sup>5</sup> ID, Steven Gillman<sup>3,6</sup> ID,
Carlos Gómez-Guijarro<sup>7</sup> ID, Katriona M. L. Gould<sup>3,4</sup>, Kasper E. Heintz<sup>3,4</sup> ID, Olivier Ilbert<sup>8</sup> ID, Christian Kragh Jespersen<sup>9</sup> ID,
Vasily Kokorev<sup>3,10</sup> ID, Mariko Kubo<sup>11</sup> ID, Georgios E. Magdis<sup>3,4,6</sup> ID, Conor J. R. McPartland<sup>3,4</sup> ID, Masato Onodera<sup>12,13</sup> ID,
Francesca Rizzo<sup>3,4</sup>, Masayuki Tanaka<sup>12,14</sup>, Sune Toft<sup>3,4</sup> ID, Aswin P. Vijayan<sup>3,6</sup>, John R. Weaver<sup>15</sup> ID,
Katherine E. Whitaker<sup>3,15</sup> ID, and Lillian Wright<sup>15,16</sup> ID

<sup>1</sup>Department of Astronomy, School of Science, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan; kei.ito@astron.s.u-tokyo.ac.jp, kei.ito.astro@gmail.com

<sup>2</sup>European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei Munchen, Germany

<sup>3</sup>Cosmic Dawn Center (DAWN), Copenhagen, Denmark

<sup>4</sup>Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark

<sup>5</sup>Caltech/IPAC, 1200 E. California Boulevard, Pasadena, CA 91125, USA

<sup>6</sup>DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark

<sup>7</sup>Université Paris-Сaclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France

<sup>8</sup>Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

<sup>9</sup>Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

<sup>10</sup>Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700AV Groningen, The Netherlands

<sup>11</sup>Astronomical Institute, Tohoku University, 6-3, Aramaki Aoba, Aoba-ku, Sendai 980-8578, Japan

<sup>12</sup>Department of Astronomical Science, The Graduate University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo, 181-8588, Japan

<sup>13</sup>Subaru Telescope, National Astronomical Observatory of Japan, National Institutes of Natural Sciences (NINS), 650 North A'ohoku Place, Hilo, HI 96720, USA

<sup>14</sup>National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo, 181-8588, Japan

<sup>15</sup>Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA

<sup>16</sup>Department of Physics and Astronomy, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul, Republic of Korea

Received 2023 July 7; revised 2024 January 31; accepted 2024 January 31; published 2024 March 29

## Abstract

We present the results of a systematic study of the rest-frame optical morphology of quiescent galaxies at $z \ge 3$ using the Near-Infrared Camera (NIRCam) on board the James Webb Space Telescope (JWST). Based on a sample selected by *UVJ* color or *NUVUVJ* color, we focus on 26 quiescent galaxies with $9.8 < \log(M_*/M_\odot) < 11.4$ at $2.8 < z_{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 and derive their size at the rest frame $0.5 \ \mu m$ using size measurements in multiple bands. Our quiescent galaxies show a significant correlation between the rest-frame $0.5 \ \mu m$ size and the stellar mass at $z \ge 3$. The analytical fit for them at $\log(M_*/M_\odot) > 10.3$ implies that our size-stellar mass relations are below those at lower redshifts, with the amplitude of $\sim 0.6 \ \text{kpc}$ at $M_* = 5 \times 10^{10} M_\odot$. This value agrees with the extrapolation of 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 mainly comprises 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 at higher redshift, suggesting $3 < z < 5$ might be the epoch of onset of morphological transformation with a fraction of very notable disky quenched galaxies.

*Unified Astronomy Thesaurus concepts: Galaxy evolution (594); Galaxy quenching (2040); High-redshift galaxies (734); Galaxy radii (617); Quenched galaxies (2016)*

## 1. Introduction

One of the most fundamental ways to characterize galaxies is to look at their appearance. 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. First, newly quenched galaxies at lower

<figure><img src="image_3.png" /></figure>

<footnote n="17"><sup>17</sup>JSPS Research Fellow (PD).</footnote>
redshift having a larger size similar to star-forming galaxies can cause an apparent size evolution, which is called progenitor bias (e.g., Barro et al. 2013; Carollo et al. 2013; Cassata et al. 2013). The dry minor mergers also effectively increase the size of quiescent galaxies (e.g., Nipoti et al. 2012; Oser et al. 2012). In addition, this trend might be only in the half-light size. Suess et al. (2019) show 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 depending on the wavelength. To minimize its effect, some studies derive the half-light radius in the same rest-frame wavelength by 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; Carnall et al. 2023a, 2023b; Nanayak-kara et al. 2024). The morphology of some of these galaxies is also explored using ground-based telescopes. 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 \gtrsim 3$ spectroscopically confirmed quiescent galaxies (Gobat et al. 2012; Saracco et al. 2020; Esdaile et al. 2021; Lustig et al. 2021).

However, a statistical discussion of the structural properties was not carried out for quiescent galaxies at $z \gtrsim 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 that resolve these galaxies at $\gtrsim 2 \mu$m, corresponding to the rest-frame wavelength of $0.5 \mu$m at $z \gtrsim 3$ often used in the size-mass relation at lower redshift. This wavelength range is not observable with HST. For $z > 3$ galaxies, HST can only probe a wavelength of $<0.4 \mu$m in the rest frame. Other approaches, such as using AO, lead to expensive observation and a severely limited and incomplete sample. Second, since the method for selecting quiescent galaxies differs among the literature above (e.g., based on the specific star formation rate or colors, also see Table 4 in Valentino et al. 2023, for a summary of the literature), comparisons between different samples are not straightforward. Third, assembling a large sample of this rare galaxy population requires large and deep imaging, especially in the 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 \gtrsim 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 the near-infrared, particularly at $\lambda > 2 \mu$m, allows us to probe the emission from the stellar populations of galaxies at $z \gtrsim 3$ and accurately measure their structural properties in the rest-frame optical. Moreover, utilizing its multiband 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 $z \sim 8$, including quiescent ones (Suess et al. 2022; Ormerod et al. 2024). The size evolution is also explored for star-forming galaxies at higher redshift (e.g., Yang et al. 2022; Gómez-Guijarro et al. 2023; Ono et al. 2023; Ward et al. 2024).

In this paper, we systematically investigate the structural properties of a sample of 26 quiescent galaxies at $z \gtrsim 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$ colors (Williams et al. 2009) or $NUV-U, U-V$, and $V-J$ colors (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. 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 \text{ km s}^{-1}\text{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 \gtrsim 3$ selected from JWST data, constructed in V23. Here, we provide a brief overview of the sample selection and refer the reader to V23 for more details.

V23 process all publicly available JWST imaging of fields spanning ~145 arcmin² taken during the first 3 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. 2024, 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 \gtrsim 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). Here, we 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),
<figure><img src="image_4.png" /><figcaption>**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 \times 4''0$ centered on each quiescent galaxy. They are sorted according to their photometric redshift. At the bottom of each figure, the ID of galaxies, identical to those of V23, and the method that selects it as quiescent one are shown. We can see that most of them have compact profiles.</figcaption></figure>

which constructs a Gaussian mixture model for the quiescent galaxy selection at high redshift (Gould 2023). Based on the *NUVUVJ* diagram, this model returns the probability of being quiescent, $P_Q$. Gould et al. (2023) show that the selection based on $P_Q$ 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$. This study selects quiescent galaxies based on the $P_{Q,50\%} > 0.7$, where $P_{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 among those in V23 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. 16 quiescent galaxies are selected from both *UVJ* and *NUVUVJ* selections, which is likely because the threshold of $P_Q$ in *NUVUVJ* selection is strict, and we select a similar population to those selected from *UVJ*.

We impose three additional selections to the sample of V23 to derive morphological properties accurately. First, 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" 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 include 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 selected only by *NUVUVJ*, leading 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 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 \le z \le 4.63$ with $9.77 \le \log(M_*/M_\odot) \le 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 spectral energy distribution (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, PI: J. S. Dunlop), 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).

### 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)
<figure><img src="image_5.png" /><figcaption>**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 UVJ and NUVUVJ. They are colored according to their photometric redshift, and the markers are larger for more massive galaxies. The gray background represents the distribution of galaxies with $log(M_*/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.</figcaption></figure>

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 in Section 3.3 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 each filter's cutouts of $4''_0 \times 4''_0$ 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, it is renormalized by a correction factor, which makes the sigma image of the regions without any sources equal to its fluctuation. Object masks are generated using the photutils module. Objects more than 1" away from the target objects are masked.

### 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 2024), 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 to the median profile of natural stars of 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. 2023) and is free from the noises that 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 (P.A.) 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 FLUXRadius 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 half-width at half-maximum (HWHM) of the used PSF is 0.''03–0.''08, which is larger at a longer wavelength (see Appendix Table 3). Appendix 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 µm, similar to van der Wel et al. (2014) for our galaxies ($\lambda_{rest} = 0.49–0.72$ µm depending on their redshift). In the fitting, the magnitude *m*, effective radius $r_{eff}$, Sérsic index *n*, axis ratio *q*, P.A., 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 \text{ mag} < m - m_{SEP} < 3 \text{ mag}$, $0.1 \text{ pixel} < r_{eff} < 500 \text{ pixel}$, $0.2 < n < 8$, and $0.0001 < q < 1$, where $m_{SEP}$ is the Kron magnitude from SEP. Next, we fit the images in the other filters. In this case, we fix the axis ratio and P.A. 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 consistently 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, but the result does not change significantly even if we fix it to the value derived independently. Also, surrounding objects within 1" from the target galaxies are simultaneously fit by the Sérsic profile.

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
<figure><img src="image_6.png" /><figcaption>**Figure 3.** A GALFIT fit of #7995 as an example, which is UVJ-selected quiescent galaxy at $z = 3.21$ with $\log(M_*/M_\odot) = 10.54$ in the CEERS field. The observed object, mask, model, and residual images 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'' \times 4.0''$. The images are scaled by the noise ($\pm 15\sigma$) of object images in each filter. The apparent size in these images does not necessarily reflect the actual size of galaxies due to the image scaling and the PSF smoothing.</figcaption></figure>

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 Sections 4 and 5 (see Section 3.5 for the selection). Some galaxies have compact residuals with both positive and negative pixels. These residuals 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 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 = 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, we conduct the following procedure. 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 mock 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 percentile values of this distribution. We find that the original uncertainty reported by GALFIT is two to five times more underestimated depending on wavelength than that derived by this procedure.

### 3.4. Bias and Quality in the Fitting

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

First, 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 \text{ mag} < m < 29 \text{ mag}, 0.002 < r_{\text{eff}} < 2''$, and $0.1 < q < 1$. The P.A. is set to be random. Since we are focusing on quiescent galaxies, which are like to have bulge-like structures as supported by the literature (e.g., Lustig et al. 2021) and Section 5.1 of this paper, this section assumes the Sérsic index as $n = 3-5$. For reference, Appendix C shows the case of $n = 1-3$. Its trend is broadly consistent with the case with $n = 3-5$, except for the smallest input size ($0.002 < r_{\text{eff}} < 0.02$). 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 even if we use images in different filters. In total, we generate 90,000 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 ~26.3 mag, which corresponds to the S/N = 50 detection in the CEERS F277W image (see V23 for the depth of images). First, the effective radii are slightly overestimated in smaller input sizes ($20\%$ in $0.002 < r_{\text{eff}} < 0.02$) and slightly underestimated in larger input sizes ($20\%$ in $0.2 < r_{\text{eff}} < 2''$). 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 < r_{\text{eff}} < 0.2$, but underestimated in larger or smaller sizes, in particular larger axis ratio ($q > 0.4$). Faintier (brighter) objects are expected to have more (less) significant scatter and bias in all parameters.

As described in the following sections, our sample typically has $r_{\text{eff}} \sim 0.6 \text{ kpc}$ ($\sim 0.01''$) and $q \sim 0.6–0.8$. Therefore, this test implies that our sample with our S/N cut of S/N > 50 is typically expected to have only a negligible offset in the effective radii and axis ratio ($<10\%$) and a slight offset in the Sérsic index ($<30\%$).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Median relative uncertainties of the estimated parameters (effective radius: left, Sérsic index: middle, axis ratio: right) derived through mock galaxy images with a Sérsic index of *n* = 3–5 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 < *r*<sub>eff</sub> < 0.02, 0.02 < *r*<sub>eff</sub> < 0.2, and 0.2 < *r*<sub>eff</sub> < 2". The shaded regions correspond to the ranges between the 16th and 84th percentile values. The vertical line in each panel corresponds to the S/N = 50 in the CEERS F277W image.</figcaption></figure>

## 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 for a fair comparison with the literature at lower redshift.

We obtain the effective radius in the rest frame 0.5 µ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_{\text{eff}}(\lambda_{\text{rest}}) = \gamma \log \left( \frac{\lambda_{\text{rest}}}{0.5 \ \mu\text{m}} \right) + \log r_{\text{eff},0.5\mu\text{m}}, \quad (1) $$

where $\gamma$ and $r_{\text{eff},0.5\mu\text{m}}$ correspond to the strength of the wavelength dependence of the effective radius ($\gamma = \Delta \log r_{\text{eff}}/\Delta \log \lambda_{\text{rest}}$) and the effective radius at the rest frame 0.5 µm. $\lambda_{\text{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_{\text{rest}}$ (~3%). Following the literature, we use the effective radius along the semimajor 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 an S/N > 50. This detection threshold is to assure the measurement quality in each band and 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_{\text{eff}} = 0.1$ pixel or 500 pixels) are not used here. After these selections, we require that at least three bands be used. The number of data points in this fitting depends on the field and the target's brightness, which spans from three to seven.

We successfully determine $\gamma$ and $r_{\text{eff},0.5\mu\text{m}}$ of 26 quiescent galaxies, which we use in the following analysis. Two galaxies (#7912 and #362) are removed from the discussion because they are too compact and have the best effective radii at the edge of the parameter range in all wavelengths at least longer than 2.7 µm and thus do not have enough size measurements. 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 samples is $\gamma = -0.38_{-0.08}^{+0.08}$ and $\gamma = -0.33_{-0.11}^{+0.10}$,
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Wavelength dependence of the effective radius ($\gamma = \Delta \log r_{\text{eff}}/\Delta \log \lambda_{\text{rest}}$) as a function of stellar mass for UVJ (left) and NUVUVJ (right) selected quiescent galaxies. #185 is not shown since its wavelength dependence is hardly constrained. Each marker is colored according to its redshift. Large squares show the median of each sample. A dashed line represents the relation of quiescent galaxies at $0.8 < z < 1.0$ in Kawinwanichakij et al. (2021).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Size–mass relation of UVJ-selected (left) and NUVUVJ-selected (right) quiescent galaxies at $z \ge 3$. The effective radii are measured at the rest frame 0.5 µm. Circles represent quiescent galaxies from this study, color-coded by their photometric redshift. The red line and shaded region represent the best fit of the size–mass relation of quiescent galaxies from this study and its 1σ uncertainty, respectively. The data points located in the gray hatched region ($\log(M_*/M_{\odot}) < 10.3$) are not used for the analytical fit. The solid and dashed black lines represent 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.</figcaption></figure>

respectively. This negative correlation between wavelength and 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.

# 4. Rest-frame Optical Size–Mass Relation

## 4.1. Overall Trends

We show quiescent galaxies at $z \ge 3$ in the size–mass plane in Figure 6. The effective radii at the rest frame 0.5 µ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 in 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_*/M_\odot) < 10.5$) have $r_{\text{eff}} > 2$ 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 and 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 (Appendix Figure 13). Also, #7399 can be affected by nearby sources. #8967 has an arclet shape, but it is not likely strongly lensed since the redshift of the companion galaxy ($z = 3.49$) is consistent with #8967. We note that removing these galaxies from the sample does not affect our conclusions. In contrast, there are galaxies with significantly small sizes ($r_{\text{eff}} < 0.1$ 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} \text{ erg cm}^{-2} \text{ s}^{-1}$ in the 0.5–2 keV. This corresponds to the limiting luminosity in the rest frame 2–10 keV of $L_X = (1.5 - 1.6) \times 10^{43} \text{ erg s}^{-1}$ at their redshift, which implies that they are at least not likely powerful unobscured X-ray AGNs. We find that removing these galaxies from the sample does not affect our conclusions, either.

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++.<fnref n="18" /> They are in good agreement with the offset of $\log(M_{*,\text{EAZY}}/M_\odot) - \log(M_{*,\text{FAST++}}/M_\odot) = 0.069^{+0.075}_{-0.066}$ (0.067$^{+0.033}_{-0.027}$) 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.

## 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; Kawinwanichakij et al. 2021; Nedkova et al. 2021). First, following van der Wel et al. (2014), the effective radii distribution is assumed to be a log-normal distribution $N(\log R_{\text{eff}}, \sigma_{\log R_{\text{eff}}})$, where $\log R_{\text{eff}}$ is the mean of the effective radius and $\sigma_{\log R_{\text{eff}}}$ is the intrinsic scatter of the distribution. We assume that the mean of the effective radius $R_{\text{eff}}$ depends on the

stellar mass as

$$R_{\text{eff}}/\text{kpc} = A \left( \frac{M_*}{5 \times 10^{10} M_\odot} \right)^\alpha, \quad (2)$$

where A is the effective radius at $M_* = 5 \times 10^{10} M_\odot$ and $\alpha$ is the slope of the size–mass relation. We consider the uncertainty $\delta r_{\text{eff}}$ of the observed effective radius $r_{\text{eff}}$ by assuming it follows a Gaussian distribution. The $\delta r_{\text{eff}}$ includes not only the observed uncertainty of $r_{\text{eff}}$ but also the uncertainty of the stellar mass, assuming 0.21 dex, the scatter of stellar mass from eazy-py (Gould et al. 2023, and also see the discussion in Section 4.1). 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_{\text{eff}}}$) as

$$\mathcal{L} = \sum \ln [W \cdot \langle N(\log R_{\text{eff}}, \sigma_{\log R_{\text{eff}}}), N(\log r_{\text{eff}}, \delta \log r_{\text{eff}}) \rangle], \quad (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. (2023) 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_*/M_\odot) = 10 - 10.5$ at $z < 2$ (e.g., Mowla et al. 2019a; Kawinwanichakij et al. 2021; Nedkova 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_*/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., ~2 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. (2023). 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_* = 5 \times 10^{10} M_\odot$), $\alpha$ (the slope of the size–mass relation), and $\sigma_{\log R_{\text{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

<footnote n="18">¹⁸ https://github.com/cschreib/fastpp</footnote>
<figure><img src="image_10.png" /><figcaption>**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 those in Figure 6. In the right panel, the best fit of the *z* < 3.3 subsample is shown as a dashed blue line for easy comparison.</figcaption></figure>

<table><caption>**Table 1**
Summary of the Obtained Best-fit Parameters of the Size–Mass Relation</caption><thead><tr><th>Sample</th><th>N<sup>a</sup></th><th>log(A/kpc)<sup>b</sup></th><th>α<sup>b</sup></th><th>σ<sub>log R<sub>eff</sub></sub><sup>b</sup></th></tr></thead><tbody><tr><td>UVJ</td><td>22</td><td>-0.18<sup>+0.06</sup><sub>-0.06</sub></td><td>0.32<sup>+0.20</sup><sub>-0.21</sub></td><td>0.18<sup>+0.07</sup><sub>-0.05</sub></td></tr><tr><td>UVJ (z < 3.3)</td><td>12</td><td>-0.14<sup>+0.08</sup><sub>-0.08</sub></td><td>0.43<sup>+0.29</sup><sub>-0.35</sub></td><td>0.16<sup>+0.09</sup><sub>-0.06</sub></td></tr><tr><td>UVJ (z > 3.3)</td><td>10</td><td>-0.26<sup>+0.15</sup><sub>-0.17</sub></td><td>0.39<sup>+0.43</sup><sub>-0.49</sub></td><td>0.28<sup>+0.13</sup><sub>-0.11</sub></td></tr><tr><td>NUVUVJ</td><td>12</td><td>-0.24<sup>+0.08</sup><sub>-0.09</sub></td><td>0.54<sup>+0.24</sup><sub>-0.28</sub></td><td>0.20<sup>+0.13</sup><sub>-0.08</sub></td></tr><tr><td>vDW14 (z = 2.75)<sup>c</sup></td><td>...</td><td>-0.06 ± 0.03</td><td>0.79 ± 0.07</td><td>0.14 ± 0.03</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The number of galaxies used in deriving the best-fit parameters of the size-mass relation.

<sup>b</sup> Their uncertainties are taken from the 16th to 84th percentiles of the total chain in the MCMC fitting.

<sup>c</sup> The best-fit parameters of early-type galaxies at 2.5 < *z* < 3.0 reported in van der Wel et al. (2014).

best-fit α confirms that the positive correlation between the stellar mass and the effective radius exists in both UVJ-selected and NUVUVJ-selected galaxies, even considering the 1σ uncertainty. The intrinsic scatter is ~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 shows 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* ≥ 3.

### 4.3. Evolution of the Size–Mass Relation of Quiescent Galaxies from *z* ~ 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 α. 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. Our measured relations are ~0.2 dex below those at *z* ~ 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}$ kpc or $R_e = 4.3 \times h(z)^{-1.29}$ kpc. 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-selected 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), particularly 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 those at lower redshift (van der Wel et al. 2014; Kawinwanichakij et al. 2021; Nedkova et al. 2021) and the size measurements of
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Left panel: redshift evolution of the effective radius at $M_{*} = 5 \times 10^{10} M_{\odot}$ in the best-fit size–mass relation. Large red and blue circles show the cases of UVJ-selected and NUVUVJ-selected quiescent galaxies. Small red circles represent the case for redshift sub-binned UVJ-selected quiescent galaxies. Black circles, gray squares, and open triangles represent 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 represent the size evolution function suggested in van der Wel et al. (2014), and the dashed lines represent their extrapolation at $z > 3$. Right panel: redshift evolution of the slope of the size–mass relation (α). The symbols’ meanings and the data points’ horizontal axis values are the same as in the left panel.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Compilation of the size–mass relation of quiescent galaxies up to $z \sim 4$. Those from this study ($z \ge 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, in van der Wel et al. (2014) in dotted lines ($0 < z < 3$), in Kawinwanichakij et al. (2021) in dashed–dotted lines ($0 < z < 1$), and in Nedkova et al. (2021) in the double-dotted line ($0 < z < 2$), respectively. Their color represents their median redshift. Markers represent the size measurement of individual quiescent galaxies from the literature at $z \ge 2.5$ (Gobat et al. 2012; Kubo et al. 2018; Saracco et al. 2020; Esdaile et al. 2021; Lustig et al. 2021; Forrest et al. 2022; Carnall et al. 2023b). They are colored by their redshift, too. We note those from Forrest et al. (2022) and Lustig et al. (2021) are corrected for the wavelength dependence of the size using the empirical relation from van der Wel et al. (2014), but the others are not. The wavelength of the used filter is different among them. If the rest-frame wavelength of the size is less than 0.4 µm, the color of their marker edge and error bar is set to be gray. It is set to be black if the rest-frame wavelength of the size is longer than 0.4 µm.</figcaption></figure>

individual quiescent galaxies at $z \ge 2.5$ in the literature (Gobat et al. 2012; Kubo et al. 2018; Saracco et al. 2020; Esdaile et al. 2021; Lustig et al. 2021; Forrest et al. 2022; Carnall et al. 2023b). Using the $K_s$-band image of the AO-assisted Very Large Telescope/High Acuity Wide-field K-band Imager (HAWK-I), we also obtain the structural parameters 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_e = 0.74_{-0.09}^{+0.20}$ kpc (see Appendix D for further details). Though the observed wavelength differs among the sample, and many of these values are not corrected for the wavelength dependence of the size (except Lustig et al. 2021; Forrest et al. 2022), they are overall consistent with our measurements. Some of the quiescent galaxies in Lustig et al. (2021) are larger than our size–mass relation, and this might be due to the redshift difference because of lower redshift for most of them than ours. It should be noted that all $z \ge 4$ measurements in Kubo et al. (2018) and Carnall et al. (2023b) are below our size–mass relations, as is SXDS-27434 at $z = 4.01$. This is in agreement with our $z \ge 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 kpc at $\log(M_{*}/M_{\odot}) = 11.0$ on average. Such a systematically small size for the $z \ge 4$ sample may imply that the size evolution of quiescent galaxies starts at $z \ge 4$.

# 5. Morphological Properties

## 5.1. Sérsic Index and Disk Fraction

The Sérsic index measured in F277W, the filter close to the rest frame 0.5 µ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., the Sérsic index is not fixed to $n = 4$ in the fitting). The measured Sérsic indices range from $n = 0.5–6$, with median values of $n = 3.42_{-0.57}^{+0.49}$ and $n = 3.61_{-0.47}^{+0.51}$ for the 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><img src="image_13.png" /><figcaption>**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; Esdaile et al. 2021; Lustig et al. 2021). The stellar mass range of their sample is log(*M*<sub>☉</sub>/*M*<sub>☉</sub>) > 10.3, log(*M*<sub>☉</sub>/*M*<sub>☉</sub>) > 11.25, 10.8 < log(*M*<sub>☉</sub>/*M*<sub>☉</sub>) < 11.3, and 11.0 < log(*M*<sub>☉</sub>/*M*<sub>☉</sub>) < 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 in van der Wel et al. (2014). All of 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*-selected and *NUVUVJ*-selected quiescent galaxies, respectively.</figcaption></figure>

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; Esdaile et al. 2021; Lustig et al. 2021). This suggests that most quiescent galaxies have bulge-dominant structures as early as *z* ~ 4.

One point that we should mention from our sample is that some quiescent galaxies have a low Sérsic index of *n* ~ 0.5–2. To investigate this further, we next derive the disk fraction, defined as the number of galaxies with *n* < 2 divided by the total number of galaxies. These values and their uncertainties 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 redshifts are consistent with zero, the disk fraction increases toward higher redshifts. In particular, the high redshift bin (z > 3.3) for *UVJ*-selected quiescent galaxies has a nonzero value (0.36<sup>+0.09</sup><sub>-0.18</sub>) even considering the

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Axis ratio in F277W filter as a function of the redshift. The small colored circles represent those of the quiescent galaxies used in this study, and the larger circles and squares represent the median of the *UVJ*-selected and *NUVUVJ*-selected samples, respectively. They are subdivided into two samples according to their redshift (z ≤ 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; Esdaile et al. 2021; Lustig et al. 2021).</figcaption></figure>

1σ 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).

## 5.2. Axis Ratio

The observed axis ratio of each galaxy projected to the images depends on its inclination angle, but we can trace the intrinsic shape of the sample by taking the median and removing the effect from the inclination angle. We find that the median axis ratio in F277W is $q = 0.61_{-0.09}^{+0.12}$ and $q = 0.78_{-0.20}^{+0.04}$ for the *UVJ*-selected and *NUVUVJ*-selected samples, respectively. These values and uncertainties 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 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 the $2.8 < z \le 3.3$ and $3.3 < z < 4.6$ samples. 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) and at $z \sim 3-4$ (Straatman et al. 2015; Esdaile et al. 2021). 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$.

# 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 V23, selecting targets from all available public JWST data taken during the first 3 months of operation 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 the observed 1.15–4.4 µm.

Thanks to the measurements in the multiple bands, we derive the rest-frame 0.5 µ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 \ge 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 the younger ones. Since a longer wavelength (i.e., $\lambda > 0.5$ µm in the rest-frame in our case) sees a more older stellar population, this can cause the observed negative 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 µ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_*/M_{\odot}) > 10.3$. These measurements suggest a correlation between the stellar mass and effective radii as early as $z \ge 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 possible redshift evolution of its amplitude at $z \ge 3$.

The effective radii at $M_* = 5 \times 10^{10} M_{\odot}$ from the best fit are ~0.6 kpc, ~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 (~0.3–0.5) than that of van der Wel et al. (2014, ~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). Several explanations exist for the possible difference between our results and van der Wel et al. (2014). First, this can be due to the increasing fraction of transitional galaxies from star forming to quiescent, particularly at lower stellar mass. For example, Suess et al. (2021) show that the median size of green valley galaxies lies between star-forming and quiescent galaxies. If this trend persists at $z \ge 3$, such larger galaxies flatten the size–mass relation. Second, 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 V23), this can be due to the small number of the sample or the cosmic variance due to our smaller survey field than that of van der Wel et al. (2014). Even though this can not fully explain the difference between our results and those of van der Wel et al. (2014) since the slope is still flatter when the sample is split into two redshift subsamples (Figure 7), this increases the size at the low stellar mass compared to high mass, due to the larger galaxies at lower

redshift, leading to the flatter slope. Lastly, the method of collecting the wavelength dependence of the size to derive the rest-frame 0.5 µm size can make the relation different. We correct it for individual galaxies as in 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 scatter (see Section 3.5), this could make the relation in van der Wel et al. (2014) different from ours. On the other hand, we have confirmed that only this can not cause this difference either by deriving the size–mass relation using the effective radius at the rest frame 0.5 µm from the observed value at F277W and the single median $\gamma$.

The high median Sérsic index ($n \sim 3–4$) and the axis ratios ($q \sim 0.6–0.8$) 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 toward higher redshift in the UVJ-selected sample. The fact that both two spectroscopically confirmed quiescent galaxies at $z > 4$ so far have low Sérsic indexes, i.e., $n = 1.23_{-0.91}^{+0.65}$ for SXDS-27434 (Appendix D) and $n = 2.3 \pm 0.3$ for those in Carnall et al. (2023b), 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 at $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 JWST multiband 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 \ge 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 \ge 4$, we can constrain the trend at an even higher redshift and obtain a statistical census into the connection between structural transformation and quenching of the current most distant quiescent galaxies.

Acknowledgments

We appreciate the anonymous referee for helpful comments and suggestions that improved the manuscript. 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 Nos. JP21K03622, JP22J00495, and JP23K13141. O.I. acknowledges the funding of the French Agence Nationale de la Recherche for the project iIMAGE (grant ANR-22-CE31-0007). The Cosmic Dawn
Center of Excellence is funded by the Danish National Research Foundation under grant DNRF140. G.E.M. acknowledges the Villum Fonden research grant 13160 “Gas to stars, stars to dust: tracing star formation across cosmic time,” grant 37440, “The Hidden Cosmos.”

Facility: 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).

# 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 an S/N > 200 in 0.5" aperture magnitude are selected to focus only 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, 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 FLUXRadius 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 radial profile of 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 more compact than the natural stars, particularly at shorter wavelengths. The PSFs smoothed by the Moffat function agree mostly well with natural stars in all filters. The median HWHM of PSFs for each filter among different fields is summarized in Table 2. The scatter among the values for different fields is 0%–11%, particularly smaller values at longer wavelengths. The effective radius at the rest frame 0.5 µm of our sample is determined primarily by the size measurement at wavelengths longer than 2 µm, due to the S/N cut for the

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Summary of PSF of F115W, F150W, F200W, F277W, F356W, F410M, and F444W filters 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'$. 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.</figcaption></figure>
<table><caption>**Table 2**
Summary of the Median HWHM of the PSFs Used for Images of Each Filter among Fields</caption><thead><tr><td>Filter Name</td><td>Wavelength (μm)</td><td>HWHM (arcsec)</td><td>Field Variation (%)</td></tr></thead><tbody><tr><td>f115w</td><td>1.15</td><td>0.033</td><td>6.6</td></tr><tr><td>f150w</td><td>1.5</td><td>0.035</td><td>11.1</td></tr><tr><td>f200w</td><td>2.0</td><td>0.039</td><td>5.1</td></tr><tr><td>f277w</td><td>2.77</td><td>0.062</td><td>4.6</td></tr><tr><td>f356w</td><td>3.56</td><td>0.069</td><td>1.5</td></tr><tr><td>f410m</td><td>4.1</td><td>0.075</td><td>1.0</td></tr><tr><td>f444w</td><td>4.44</td><td>0.079</td><td>0.2</td></tr></tbody></table>

magnitude. At that wavelength range, the PSFs are determined homogeneously among fields.

**Appendix B**
**Summary of Fitting Result**

Figures 13–15 show the fitting results of all 26 objects used in this study. F200W and F277W are selected here because they correspond to the rest frame ~0.5 µm at the redshift of our interest. In addition, Table 3 summarizes the best-fit parameters obtained in fitting F277W images, the effective radius at rest-frame 0.5 µm, and the wavelength dependence of the size derived in Section 3.5.
<figure><img src="image_16.png" /><figcaption>**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 \times 4''0$. The images are scaled by their noise (±15σ). The black contours start at 10% of the peak pixel count of each object and increase by $\sqrt{2}$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 14. Continued from Figure 13.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 15. Continued from Figure 13.</figcaption></figure>
<table><caption>Table 3
Summary of GALFIT Fitting</caption><thead><tr><th>ID<sup>a</sup></th><th>Field</th><th>z<sub>phot</sub><sup>b</sup></th><th>log(M<sub>*</sub>/M<sub>☉</sub>)<sup>b</sup></th><th>r<sub>eff,f277w</sub><br>(kpc, arcsec)</th><th>n<sub>f277w</sub></th><th>q<sub>f277w</sub></th><th>P.A.<sub>f277w</sub><br>(deg)</th><th>m<sub>f277w</sub><br>(mag)</th><th>r<sub>eff,0.5µm</sub><sup>c</sup><br>(kpc)</th><th>γ<sup>d</sup></th><th>N<sub>filter</sub><sup>e</sup></th><th>flag<sup>f</sup></th></tr></thead><tbody><tr><td>7174</td><td>CEERS</td><td>2.83</td><td>10.34</td><td>0.621 ± 0.021 (0.079 ± 0.003)</td><td>2.25 ± 0.23</td><td>0.701 ± 0.025</td><td>-61.0 ± 3.0</td><td>24.25 ± 0.02</td><td>0.798 ± 0.057</td><td>-0.57 ± 0.11</td><td>5</td><td>0</td></tr><tr><td>19471</td><td>J1235</td><td>2.83</td><td>10.46</td><td>0.269 ± 0.006 (0.034±0.001<sup>g</sup>)</td><td>5.65 ± 0.4</td><td>0.315 ± 0.012</td><td>14.0 ± 1.0</td><td>23.44 ± 0.01</td><td>0.464 ± 0.064</td><td>-0.87 ± 0.3</td><td>5</td><td>0</td></tr><tr><td>788</td><td>SGAS1723</td><td>2.85</td><td>10.72</td><td>3.744 ± 0.988 (0.479 ± 0.126)</td><td>4.57 ± 0.96</td><td>0.365 ± 0.034</td><td>31.0 ± 3.0</td><td>22.55 ± 0.16</td><td>4.407 ± 3.339</td><td>-0.73 ± 1.02</td><td>3</td><td>0</td></tr><tr><td>2278</td><td>PRIMER</td><td>2.85</td><td>10.41</td><td>0.722 ± 0.094 (0.092 ± 0.012)</td><td>5.52 ± 0.97</td><td>0.844 ± 0.055</td><td>-78.0 ± 15.0</td><td>23.76 ± 0.06</td><td>0.893 ± 0.279</td><td>-0.49 ± 0.44</td><td>4</td><td>0</td></tr><tr><td>4675</td><td>CEERS</td><td>2.91</td><td>10.78</td><td>0.5 ± 0.012 (0.064 ± 0.002)</td><td>4.78 ± 0.35</td><td>0.514 ± 0.01</td><td>73.0 ± 1.0</td><td>22.82 ± 0.01</td><td>0.602 ± 0.032</td><td>-0.51 ± 0.08</td><td>7</td><td>0</td></tr><tr><td>2809</td><td>PRIMER</td><td>2.92</td><td>10.77</td><td>0.823 ± 0.035 (0.106 ± 0.005)</td><td>2.89 ± 0.27</td><td>0.508 ± 0.016</td><td>77.0 ± 1.0</td><td>22.96 ± 0.02</td><td>0.926 ± 0.066</td><td>-0.3 ± 0.1</td><td>6</td><td>0</td></tr><tr><td>28463</td><td>QUINTET</td><td>2.93</td><td>10.4</td><td>1.023 ± 0.2 (0.132 ± 0.026)</td><td>[4.0]</td><td>0.626 ± 0.08</td><td>-32.0 ± 8.0</td><td>24.49 ± 0.12</td><td>0.806 ± 0.402</td><td>0.82 ± 0.74</td><td>3</td><td>2</td></tr><tr><td>2560</td><td>PRIMER</td><td>2.93</td><td>11.24</td><td>1.587 ± 0.037 (0.205 ± 0.005)</td><td>2.25 ± 0.08</td><td>0.843 ± 0.01</td><td>-45.0 ± 2.0</td><td>21.78 ± 0.01</td><td>1.836 ± 0.113</td><td>-0.34 ± 0.09</td><td>6</td><td>0</td></tr><tr><td>7399</td><td>CEERS</td><td>2.98</td><td>10.21</td><td>2.508 ± 0.202 (0.325 ± 0.026)</td><td>[4.0]</td><td>0.243 ± 0.02</td><td>37.0 ± 1.0</td><td>24.5 ± 0.06</td><td>4.882 ± 1.185</td><td>-1.74 ± 0.38</td><td>5</td><td>2</td></tr><tr><td>277</td><td>PRIMER</td><td>3.12</td><td>10.62</td><td>0.336 ± 0.012 (0.044±0.002<sup>g</sup>)</td><td>2.11 ± 0.31</td><td>0.86 ± 0.043</td><td>35.0 ± 11.0</td><td>23.39 ± 0.01</td><td>0.35 ± 0.021</td><td>-0.07 ± 0.11</td><td>5</td><td>0</td></tr><tr><td>2796</td><td>JWST-NEP-TDF-NRC</td><td>3.13</td><td>10.53</td><td>0.809 ± 0.047 (0.106 ± 0.006)</td><td>3.88 ± 0.35</td><td>0.817 ± 0.02</td><td>75.0 ± 3.0</td><td>23.38 ± 0.03</td><td>0.952 ± 0.103</td><td>-0.65 ± 0.17</td><td>5</td><td>0</td></tr><tr><td>11266</td><td>CEERS</td><td>3.18</td><td>10.74</td><td>0.494 ± 0.018 (0.065 ± 0.002)</td><td>4.89 ± 0.41</td><td>0.862 ± 0.017</td><td>61.0 ± 4.0</td><td>23.12 ± 0.02</td><td>0.575 ± 0.046</td><td>-0.44 ± 0.13</td><td>6</td><td>0</td></tr><tr><td>7955</td><td>CEERS</td><td>3.21</td><td>10.54</td><td>0.517 ± 0.011 (0.069 ± 0.001)</td><td>3.52 ± 0.26</td><td>0.566 ± 0.012</td><td>77.0 ± 1.0</td><td>23.14 ± 0.01</td><td>0.591 ± 0.046</td><td>-0.57 ± 0.14</td><td>6</td><td>0</td></tr><tr><td>9621</td><td>CEERS</td><td>3.38</td><td>10.98</td><td>1.272 ± 0.015 (0.172 ± 0.002)</td><td>1.7 ± 0.05</td><td>0.421 ± 0.006</td><td>9.0 ± 0.0</td><td>22.9 ± 0.01</td><td>1.374 ± 0.028</td><td>-0.37 ± 0.04</td><td>6</td><td>0</td></tr><tr><td>15782</td><td>CEERS</td><td>3.35</td><td>11.39</td><td>1.711 ± 0.043 (0.23 ± 0.006)</td><td>4.05 ± 0.1</td><td>0.806 ± 0.005</td><td>79.0 ± 1.0</td><td>21.34 ± 0.01</td><td>1.857 ± 0.073</td><td>-0.38 ± 0.07</td><td>7</td><td>0</td></tr><tr><td>35943</td><td>QUINTET</td><td>3.36</td><td>10.12</td><td>0.708 ± 0.067 (0.095 ± 0.009)</td><td>1.98 ± 0.43</td><td>0.905 ± 0.063</td><td>-81.0 ± 82.0</td><td>24.5 ± 0.05</td><td>0.719 ± 0.134</td><td>-0.11 ± 0.32</td><td>3</td><td>0</td></tr><tr><td>8083</td><td>CEERS</td><td>3.41</td><td>10.87</td><td>0.672 ± 0.013 (0.091 ± 0.002)</td><td>3.17 ± 0.13</td><td>0.775 ± 0.009</td><td>-71.0 ± 1.0</td><td>22.49 ± 0.01</td><td>0.719 ± 0.023</td><td>-0.31 ± 0.06</td><td>7</td><td>0</td></tr><tr><td>1863</td><td>PRIMER</td><td>3.46</td><td>10.72</td><td>0.739 ± 0.043 (0.101 ± 0.006)</td><td>2.58 ± 0.38</td><td>0.447 ± 0.029</td><td>-28.0 ± 2.0</td><td>23.81 ± 0.04</td><td>0.848 ± 0.126</td><td>-0.41 ± 0.26</td><td>4</td><td>1</td></tr><tr><td>8967</td><td>JWST-NEP-TDF-NRC</td><td>3.47</td><td>9.77</td><td>2.214 ± 0.153 (0.301 ± 0.021)</td><td>0.52 ± 0.1</td><td>0.502 ± 0.035</td><td>60.0 ± 3.0</td><td>25.46 ± 0.09</td><td>2.159 ± 0.635</td><td>-0.37 ± 0.56</td><td>3</td><td>0</td></tr><tr><td>8782</td><td>CEERS</td><td>3.53</td><td>10.39</td><td>0.316 ± 0.012 (0.043±0.002<sup>g</sup>)</td><td>4.85 ± 0.53</td><td>0.941 ± 0.028</td><td>46.0 ± 16.0</td><td>23.62 ± 0.02</td><td>0.332 ± 0.031</td><td>-0.18 ± 0.18</td><td>6</td><td>0</td></tr><tr><td>9622</td><td>CEERS</td><td>3.54</td><td>11.11</td><td>0.487 ± 0.005 (0.067 ± 0.001)</td><td>3.64 ± 0.14</td><td>0.485 ± 0.006</td><td>60.0 ± 0.0</td><td>22.16 ± 0.01</td><td>0.532 ± 0.031</td><td>-0.11 ± 0.12</td><td>7</td><td>0</td></tr><tr><td>1427</td><td>PRIMER</td><td>3.55</td><td>10.18</td><td>0.086 ± 0.016 (0.012 ± 0.002)</td><td>[4.0]</td><td>0.514 ± 0.205</td><td>-37.0 ± 13.0</td><td>24.6 ± 0.03</td><td>0.054 ± 0.024</td><td>0.67 ± 0.95</td><td>5</td><td>3</td></tr><tr><td>10378</td><td>CEERS</td><td>3.67</td><td>10.83</td><td>1.002 ± 0.081 (0.139 ± 0.011)</td><td>5.34 ± 0.53</td><td>0.648 ± 0.016</td><td>70.0 ± 2.0</td><td>23.26 ± 0.04</td><td>1.161 ± 0.106</td><td>-0.85 ± 0.17</td><td>6</td><td>0</td></tr><tr><td>2876</td><td>CEERS</td><td>4.15</td><td>10.81</td><td>0.722 ± 0.014 (0.105 ± 0.002)</td><td>1.67 ± 0.11</td><td>0.464 ± 0.016</td><td>17.0 ± 1.0</td><td>23.98 ± 0.01</td><td>0.743 ± 0.019</td><td>-0.45 ± 0.06</td><td>5</td><td>0</td></tr><tr><td>185</td><td>PRIMER</td><td>4.51</td><td>10.91</td><td>0.033 ± 0.035 (0.005 ± 0.005)</td><td>4.66 ± 3.52</td><td>0.232 ± 0.313</td><td>-82.0 ± 81.0</td><td>23.94 ± 0.09</td><td>0.034 ± 0.046</td><td>0±4 × 10<sup>7</sup></td><td>3</td><td>1</td></tr><tr><td>10084</td><td>CEERS</td><td>4.63</td><td>11.23</td><td>0.579 ± 0.016 (0.089 ± 0.002)</td><td>[4.0]</td><td>0.968 ± 0.015</td><td>34.0 ± 16.0</td><td>22.95 ± 0.01</td><td>0.475 ± 0.038</td><td>0.15 ± 0.21</td><td>7</td><td>2</td></tr></tbody></table>

Notes.

a The IDs are identical to those in V23.

b They are derived with eazy-py in V23.

c The effective radius at the rest frame 0.5 µm derived in Section 3.5.

d The wavelength dependence of the effective radius derived in Section 3.5.

e Number of filters used in deriving the γ and r<sub>eff,0.5µm</sub>.

f The quality of fitting. flag =0: good fit. flag =1: 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.

g Though the measured effective radii are smaller than the HWHM of PSFs in F277W, they have larger effective radii than the HWHM of PSFs in the other filters at shorter wavelengths.
## Appendix C
### Bias and the Quality in the Fitting for the Profile with $n = 1-3$

Figure 16 shows the relative uncertainties of the effective radius, Sérsic index, and axis ratio for a mock Sérsic profile

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Median relative uncertainties of the estimated parameters (effective radius: left, Sérsic index: middle, axis ratio: right) derived through mock galaxy images with Sérsic index of $n = 1-3$ as a function of the magnitude. The meanings of each panel and the lines are the same as in Figure 4.</figcaption></figure>

with $n = 1-3$. The trend of the profile with the effective radius of $r_e > 0''02$ is overall consistent with those with $n = 3-5$. It is not the case for $r_e < 0''02$, where the effective radius tends to be slightly smaller and the Sérsic index tends to be larger than the input values.

## Appendix D
### New Size Estimate of SXDS-27434 at $z = 4.01$

Here, we report the size estimate at rest frame ~4200 Å of the galaxy SXDS-27434 at $z = 4.01$ (Tanaka et al. 2019; Valentino et al. 2020). This is based on observations with 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 (R. A., decl. = 34.29871, -4.98987) was observed with the Ks filter for 2040 s on source in service mode on 2019 October 22. We applied the standard jitter scheme within a 25" box and obtained 34 offset exposures (NDIT × DIT = 6 × 10 s). We reduce the data as described in Brammer et al. (2016).<fnref n="19" /> 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''08, and their flux calibration is 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 an FWHM = 0''34. The best-fit model and the residuals are shown in Figure 17. We measure an effective semimajor axis of $r_e = 0.74_{-0.09}^{+0.20}$ kpc and a Sérsic index $n = 1.23_{-0.91}^{+0.65}$ with an axis ratio $q = 0.88_{-0.26}^{+0.06}$. The uncertainties are estimated by bootstrapping 1000× 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 16th to 84th 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''0.8 or 0''1) and a fixed or freely varying sky level. The target is resolved, and the

<footnote n="19">¹⁹ https://github.com/gbrammer/HAWKI-FF</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure 17.** 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.</figcaption></figure>

radius is 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_e = 0.76 \pm 0.20$ kpc).

## ORCID iDs

Kei Ito https://orcid.org/0000-0002-9453-0381

Francesco Valentino https://orcid.org/0000-0001-6477-4011

Gabriel Brammer https://orcid.org/0000-0003-2680-005X

Andreas L. Faisst https://orcid.org/0000-0002-9382-9832

Steven Gillman https://orcid.org/0000-0001-9885-4589

Carlos Gómez-Guijarro https://orcid.org/0000-0002-4085-9165

Kasper E. Heintz https://orcid.org/0000-0002-9389-7413

Olivier Ilbert https://orcid.org/0000-0002-7303-4397

Christian Kragh Jespersen https://orcid.org/0000-0002-8896-6496

Vasily Kokorev https://orcid.org/0000-0002-5588-9156

Mariko Kubo https://orcid.org/0000-0002-7598-5292

Georgios E. Magdis https://orcid.org/0000-0002-4872-2294

Conor J. R. McPartland https://orcid.org/0000-0003-0639-025X

Masato Onodera https://orcid.org/0000-0003-3228-7264

Sune Toft https://orcid.org/0000-0003-3631-7176

John R. Weaver https://orcid.org/0000-0003-1614-196X

Katherine E. Whitaker https://orcid.org/0000-0001-7160-3632

Lillian Wright https://orcid.org/0000-0002-3722-322X

## References

Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJL, 946, L12

Barbary, K. 2016, JOSS, 1, 58

Barro, G., Faber, S. M., Pérez-González, P. G., et al. 2013, ApJ, 765, 104

Bertin, E. 2011, in ASP Conf. Ser. 442, Astronomical Data Analysis Software and Systems XX, ed. I. N. Evans et al. (San Francisco, CA: ASP), 435

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bezanson, R., van Dokkum, P. G., Tal, T., et al. 2009, ApJ, 697, 1290

Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, astropy/photutils: 1.5.0, v1.5.0, Zenodo, doi:10.5281/zenodo.6825092

Brammer, G. 2021, eazy-py, v0.5.2, Zenodo, doi:10.5281/zenodo.7575984

Brammer, G., & Matharu, J. 2021, gbrammer/grizli: Release 2021, v1.3.2, Zenodo, doi:10.5281/zenodo.5012699

Brammer, G., Strait, V., Matharu, J., & Momcheva, I. 2022, grizli, v1.5.0, Zenodo, doi:10.5281/zenodo.6672538

Brammer, G. B., Marchesini, D., Labbé, I., et al. 2016, ApJS, 226, 6

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Carnall, A. C., McLeod, D. J., McLure, R. J., et al. 2023a, MNRAS, 520, 3974

Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2023b, Natur, 619, 716

Carollo, C. M., Bschorr, T. J., Renzini, A., et al. 2013, ApJ, 773, 112

Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2022, ApJ, 954, 31

Cassata, P., Giavalisco, M., Williams, C. C., et al. 2013, ApJ, 775, 106

Chabrier, G. 2003, PASP, 115, 763

Cutler, S. E., Whitaker, K. E., Mowla, L. A., et al. 2022, ApJ, 925, 34

D'Eugenio, C., Daddi, E., Gobat, R., et al. 2021, A&A, 653, A32

Ding, X., Silverman, J. D., & Onoue, M. 2022, ApJL, 939, L28

Esdaile, J., Glazbrook, K., Labbé, I., et al. 2021, ApJL, 908, L35

Faisst, A. L., Carollo, C. M., Capak, P. L., et al. 2017, ApJ, 839, 71

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Forrest, B., Marsan, Z. C., Annunziatella, M., et al. 2020, ApJ, 903, 47

Forrest, B., Wilson, G., Muzzin, A., et al. 2022, ApJ, 938, 109

Fruchter, A. S., & Hook, R. N. 2002, PASP, 114, 144

Fudamoto, Y., Inoue, A. K., & Sugahara, Y. 2022, ApJL, 938, L24

Gillman, S., Gullberg, B., Brammer, G., et al. 2023, A&A, 676, A26

Glazbrook, K., Schreiber, C., Labbé, I., et al. 2017, Natur, 544, 71

Gobat, R., Strazzullo, V., Daddi, E., et al. 2012, ApJL, 759, L44

Gómez-Guijarro, C., Magnelli, B., Elbaz, D., et al. 2023, A&A, 677, A34

Gould, K. M. L. 2023, GMM-quiescent, v1.0.0, Zenodo, doi:10.5281/zenodo.7712937

Gould, K. M. L., Brammer, G., Valentino, F., et al. 2023, AJ, 165, 248

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hunter, J. D. 2007, CSE, 9, 90

Kawinwanichakij, L., Silverman, J. D., Ding, X., et al. 2021, ApJ, 921, 38

Kissler-Patig, M., Pirard, J. F., Casali, M., et al. 2008, A&A, 491, 941

Kocevski, D. D., Hasinger, G., Brightman, M., et al. 2018, ApJS, 236, 48

Kubo, M., Tanaka, M., Yabe, K., et al. 2018, ApJ, 867, 1

Lam, S. K., Pitrou, A., & Seibert, S. 2015, in Proc. Second Workshop on the LLVM Compiler Infrastructure in HPC - LLVM '15 (New York: ACM), 1

Lustig, P., Strazzullo, V., D'Eugenio, C., et al. 2021, MNRAS, 501, 2659

Marchesini, D., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2009, ApJ, 701, 1765

Marsan, Z. C., Marchesini, D., Muzzin, A., et al. 2019, ApJ, 871, 201

McKinney, W. 2010, in Proc. IXth Python in Science Conf. (SciPy 2010) (Austin, TX: SciPy)

Mosleh, M., Tacchella, S., Renzini, A., et al. 2017, ApJ, 837, 2

Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113

Mowla, L., van der Wel, A., van Dokkum, P., & Miller, T. B. 2019a, ApJL, 872, L13

Mowla, L. A., van Dokkum, P., Brammer, G. B., et al. 2019b, ApJ, 880, 57

Nanayakkara, T., Glazbrook, K., Jacobs, C., et al. 2024, NatSR, 14, 3724

Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, ApJS, 220, 10

Nedkova, K. V., Häußler, B., Marchesini, D., et al. 2021, MNRAS, 506, 928

Newman, A. B., Belli, S., Ellis, R. S., & Patel, S. G. 2018, ApJ, 862, 126
Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, v0.8.0, Zenodo, doi:10.5281/zenodo.11813

Nipoti, C., Treu, T., Leauthaud, A., et al. 2012, MNRAS, 422, 1714

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Ono, Y., Harikane, Y., Ouchi, M., et al. 2023, ApJ, 951, 72

Ormerod, K., Conselice, C. J., Adams, N. J., et al. 2024, MNRAS, 527, 6110

Oser, L., Naab, T., Ostriker, J. P., & Johansson, P. H. 2012, ApJ, 744, 63

Pacifici, C., Iyer, K. G., Mobasher, B., et al. 2023, ApJ, 944, 141

Paufique, J., Bruton, A., Glindemann, A., et al. 2010, Proc. SPIE, 7736, 77361P

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097

Pérez-González, P. G., Barro, G., Annunziatella, M., et al. 2023, ApJL, 946, L16

Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, Proc. SPIE, 9143, 91433X

Perrin, M. D., Soummer, R., Elliott, E. M., Lallo, M. D., & Sivaramakrishnan, A. 2012, Proc. SPIE, 8442, 84423D

Pirard, J.-F., Kissler-Patig, M., Moorwood, A., et al. 2004, Proc. SPIE, 5492, 1763

Pontoppidan, K. M., Barrientes, J., Blome, C., et al. 2022, ApJL, 936, L14

Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123

Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, A&A, 558, A33

Saracco, P., Marchesini, D., La Barbera, F., et al. 2020, ApJ, 905, 40

Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018, A&A, 618, A85

Sérsic, J. L. 1963, BAAA, 6, 41

Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978

Steinhardt, C. L., Jespersen, C. K., & Linzer, N. B. 2021, ApJ, 923, 8

Straatman, C. M. S., Labbé, I., Spitler, L. R., et al. 2015, ApJL, 808, L29

Suess, K. A., Bezanson, R., Nelson, E. J., et al. 2022, ApJL, 937, L33

Suess, K. A., Kriek, M., Price, S. H., & Barro, G. 2019, ApJ, 877, 103

Suess, K. A., Kriek, M., Price, S. H., & Barro, G. 2021, ApJ, 915, 87

Tanaka, M., Valentino, F., Toft, S., et al. 2019, ApJL, 885, L34

Toft, S., Zabl, J., Richard, J., et al. 2017, Natur, 546, 510

Tzanavaris, P., Gallagher, S. C., Hornschemeier, A. E., et al. 2014, ApJS, 212, 9

Valentino, F. 2023, "An atlas of color-selected quiescent galaxies at z>3 in public JWST fields" (Supplementary material), v1.0, Zenodo, doi:10.5281/zenodo.7614908

Valentino, F., Brammer, G., Gould, K. M. L., et al. 2023, ApJ, 947, 20

Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93

van der Wel, A., Bell, E. F., Häussler, B., et al. 2012, ApJS, 203, 24

van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28

van Dokkum, P. G., Nelson, E. J., Franx, M., et al. 2015, ApJ, 813, 23

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Ward, E. M., de la Vega, A., Mobasher, B., et al. 2024, ApJ, 962, 176

Weaver, J. R., Davidzon, I., Toft, S., et al. 2023, A&A, 677, A184

Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879

Yang, L., Morishita, T., Leethochawalit, N., et al. 2022, ApJL, 938, L17

Zhuang, M.-Y., & Shen, Y. 2024, ApJ, 962, 139