Size–Stellar Mass Relation and Morphology of Quiescent Galaxies at z ≥ 3 in Public JWST Fields

We present the results of a systematic study of the rest-frame optical morphology of quiescent galaxies at z ≥ 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⊙)<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 μm using size measurements in multiple bands. Our quiescent galaxies show a significant correlation between the rest-frame 0.5 μm size and the stellar mass at z ≥ 3. The analytical fit for them at log(M⋆/M⊙)>10.3 implies that our size–stellar mass relations are below those at lower redshifts, with the amplitude of ∼0.6 kpc at M ⋆ = 5 × 1010 M ⊙. 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 ∼ 3–5. Our sample mainly comprises galaxies with bulge-like structures according to their median Sérsic index and axis ratio of n ∼ 3–4 and q ∼ 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.


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 ∼ 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 ∝ (1 + z) −1.48 at z < 3, whereas that of star-forming galaxies evolves as ∝ (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 depending 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 ∼ 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 * JSPS Research Fellow (PD) 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.
However, a statistical discussion of the structural properties was not carried out for quiescent galaxies at z ≥ 3, including the constraint on their size-mass relation and its redshift evolution to z ∼ 0. There are four underlying reasons.One is the requirement for high-resolution images which resolve these galaxies at ≥ 2 µm, corresponding to the rest-frame wavelength of 0.5 µm at z ≥ 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 the wavelength of < 0.4 µ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 literature above (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 ≥ 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 (JW ST ) is poised to revolutionize this topic.JW ST 's high-resolution imaging at near-infrared, particularly at λ > 2 µm, al-lows us to probe the emission from the stellar populations of galaxies at z ≥ 3 and accurately measure their structural 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 JW ST , such as for a quiescent galaxy at z = 4.658 (Carnall et al. 2023b), and a statistical sample of galaxies at to z ∼ 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., Ono et al. 2022;Yang et al. 2022;Gómez-Guijarro et al. 2023;Ward et al. 2023).
In this paper, we systematically investigate the structural properties of a sample of 26 quiescent galaxies at z ≥ 3 from Valentino et al. (2023) (hereafter denoted as V23) by using the JW ST /NIRCam images.V23 select quiescent galaxies in fields with publicly available JW ST data based on U − V and V − J color (Williams et al. 2009) or N U V −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.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 ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7.The magnitude is based on the AB magnitude system (Oke & Gunn 1983).We assume the initial mass function of Chabrier (2003).

SAMPLE SELECTION
This study uses the sample of quiescent galaxies at z ≥ 3 selected from JW ST data, constructed in V23.Here, we briefly overview the sample selection and refer the reader to V23 for more details.
V23 process the all publicly available JW ST imaging of fields spanning ∼ 145 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 JW ST 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 σ N M AD = 0.0268 in δz/(1+z) for the spectroscopic sample.We select galaxies at z ≥ 3 from their photometric redshift.
Two color selections for quiescent galaxies are employed.One is based on the traditional U V J 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 N U V − U , U − V , and V − J color diagram, hereafter called the N U V U V J 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).Based on the N U V U V J diagram, this model returns the "probability of being quiescent", P Q .Gould et al. (2023) shows that the selection based on P Q not only can reduce the contamination compared to the U V J 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 U V J and N U V U V J colors are the most conservative ones among those in V23 not to be affected by the contaminants due to the uncertainty in color.These U V J and N U V U V J color selections give 34 and 18 quiescent galaxies, respectively.16 quiescent galaxies are selected from both U V J and N U V U V J selections, which is likely because the threshold of P Q in N U V U V J selection is strict, and we select a similar population to those selected from U V J.
We impose three additional selections to the sample of V23 to derive morphological properties accurately.Firstly, we require galaxies to be detected in more than three bands of JW ST 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 U V J or N U V U V J, 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 N U V U V J, leading to 30 and 16 galaxies selected by U V J or N U V U V J, 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 U V J selection or N U V U V J selection.All of these galaxies are selected by U V J, and 14 galaxies are selected by N U V U V J.They are at 2.83 ≤ z ≤ 4.63 with 9.77 ≤ log (M ⋆ /M ⊙ ) ≤ 11.39 according to eazypy.Figure 1 summarizes the RGB images of the long wavelength of NIRCam for these galaxies.Their location in the U V J and N U V U V J 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)

MORPHOLOGICAL FITTING
To obtain the structural parameters, we fit an analytical function to the JW ST /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 in Section 3.3 implies that this assumption is reasonable.

Input Images
We use all available images of JW ST /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 ′′ ×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.

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(Perrin et al. , 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 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. 2022) 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.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 medianstacked natural stars.The half width at half maximum of the used PSF is 0.03 ′′ − 0.08 ′′ , which is larger at a longer wavelength (see Table 3).Figure 12 shows the one-dimensional radial profiles of smoothed PSFs used in this study in CEERS as an example.

Fitting
In fitting the Sérsic profile to images, we use Galfit (Peng et al. 2002(Peng et al. , 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 (λ 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, 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 mag < m − m SEP < 3 mag, 0.1 pixel < r eff < 500 pixel, 0.2 < n < 8, and 0.0001 < q < 1, where m SEP is the Kron magnitude from the sep.Next, we fit the images in the other filters.In this case, we fix the axis ratio and position angle as the best-fit value in the 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 ′′ × 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.(Williams et al. 2009).Squares represent quiescent galaxies selected by the U V J only, whereas circles represent those selected by U V J and N U V U V J.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⊙) > 9.5 at 3 < z < 5. Right panel: N U V U V J diagram (Gould et al. 2023).The meanings of the markers are the same as in the left panel.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 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.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 χ 2 ν = 1.0− 1.2 in each filter, suggesting that the assumption of the single Sérsic profile is overall good.

F277W images as in
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 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.

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 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.
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 mag < m < 29 mag, 0.002 ′′ < r eff < 2 ′′ , and 0.1 < q < 1.The position angle is set to be random.Since we are focusing on quiescent galaxies, which are like to have bulge-like structures supported by the literature (e.g., Lustig et al. 2021) and Section 5.1 of this paper, this section assumes 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 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 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 ∼ 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 ′′ < r eff < 0.02 ′′ ) and slightly underestimated in larger input sizes (20% in 0.2 ′′ < r eff < 2 ′′ ).The Sérsic index is underestimated (10 − 40%, equal to δn ∼ 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 eff < 0.2 ′′ , 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 typically has r eff ∼ 0.6 kpc (∼ 0.1 ′′ ) and q ∼ 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%).

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: where γ and r eff,0.5µmcorrespond to the strength of the wavelength dependence of the effective radius (γ = ∆ log r eff /∆ log λ rest ) and the effective radius at the restframe 0.5 µm.λ 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 λ rest (∼ 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 to assure the measurement quality in each band and 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 eff = 0.1 pixel or 500 pixel) 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 γ and r eff,0.5µ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 number of size measurements.Among 26 galaxies, all of them are selected by the U V J diagram, and 13 of them are selected by the N U V U V J 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 γ 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 U V J-selected and N U V U V J-selected sample is γ = −0.38 +0.08 −0.08 and γ = −0.33 +0.10 −0.11 , 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 γ = −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.

Overall Trends
We show quiescent galaxies at z ≥ 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 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 = 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 eff < 0.02 ′′ , 0.02 ′′ < r eff < 0.2 ′′ , and 0.2 ′′ < r eff < 2 ′′ .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.
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 ⊙ ) < 10.5) have r eff > 2 kpc in the U V J-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 (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 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 × 10 −16 erg/cm 2 /s 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) × 10 43 erg/s 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   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.

Analytic fits
Motivated by the appearance of the existence of sizemass 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 eff , σ log R eff ), where log R eff is the mean of the effective radius and σ log R eff is the intrinsic scatter of the distribution.We assume that the mean of the effective radius R eff depends on the stellar mass as: where A is the effective radius at M ⋆ = 5 × 10 10 M ⊙ and α is the slope of the size-mass relation.We consider the uncertainty δr eff of the observed effective radius r eff by assuming it follows a Gaussian distribution.The δr eff includes not only the observed uncertainty of r eff but also (3) where W is the weight to avoid the dominance of lowermass 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 ⋆ /M ⊙ ) = 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 ⋆ /M ⊙ ) > 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 ∼ 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 U V J 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σ 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 bestfit values of A (the effective radius at M ⋆ = 5×10 10 M ⊙ ), α (the slope of the size-mass relation), and σ log R 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 α confirms that the positive correlation between the stellar mass and the effective radius exists in both U V J-selected and N U V U V J-selected galaxies, even considering the 1σ uncertainty.The intrinsic scatter is ∼ 0.2 dex for both the U V J-selected and N U V U V J-selected sample, similar to those reported at z < 3 (e.g., van der Wel et al. 2014;Kawinwanichakij et al. 2021).
We split U V J-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 U V J-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.

Evolution of the size-mass relation of quiescent
galaxies from z ∼ 4 Figure 8 shows the evolution of the effective radius at M ⋆ = 5 × 10 10 M ⊙ (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 U V J diagram.We confirm that the effective radius at M ⋆ = 5 × 10 10 M ⊙ 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 × (1 + z) −1.48 kpc or R e = 4.3 × h(z) −1.29 kpc.The U V J-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 N U V U V Jand U V J-selected quiescent galaxiesalso 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 U V J-selected sample.The possible origin of this difference will be discussed in Section 6.
Lastly, Figure 9 compares our size-mass relation with the those at lower redshift (van der Wel et al. 2014;Kawinwanichakij et al. 2021;Nedkova et al. 2021) and the size measurements of individual quiescent galaxies at z ≥ 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 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.20  −0.09 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 ≥ 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 ≥ 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 ⊙ ) = 11.0 on average.Such a systematically small size for the z ≥ 4 sample may imply that the size evolution of quiescent galaxies starts at z ≥ 4.

Sérsic Index and Disk Fraction
The Sérsic index measured in F277W, the filter close to rest-frame 0.5 µm at z ∼ 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., Sérsic index is not fixed to n = 4 in the fitting).The measured Sérsic indices range from n = 0.5 to n = 6, with median values of n = 3.42 +0.49−0.57and n = 3.61 +0.51  −0.47 for U V J-selected and N U V U V J-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 ∼ 4.
One point that we should mention from our sample is that some quiescent galaxies have a low Sérsic index as 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 towards higher redshifts.In particular, the high redshift bin (z > 3.3) for U V J-selected quiescent galaxies has a non-zero value (0.36 +0.09 −0.18 ) even considering the 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).

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.12 −0.09 and q = 0.78 +0.04  −0.20 for U V J-selected and N U V U V Jselected sample, respectively.These values and uncertainties are again obtained from bootstrapping, as for the median value of the Sérsic index (Section 5.1).The N U V U V J-selected sample has a higher median value than the U V J-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 ≤ 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) and at z ∼ 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 ∼ 4.

DISCUSSION AND SUMMARY
This study investigates the morphological properties of 26 quiescent galaxies at 3 < z < 5 in public JW ST fields.Quiescent galaxies are identified by Valentino et al. (2023), selecting targets from all available public JW ST data taken during the first three months of operation based on two rest-frame color selections, i.e., the traditional U V J diagram (Williams et al. 2009) and a novel selection from the rest-frame N U V , 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 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 ≥ 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 longer wavelength (i.e., λ > 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 ⊙ ) > 10.3.These measurements suggest a correlation between the stellar mass and effective radii as early as z ≥ 3 with the power law slope of 0.3−0.5.In addition, by splitting the U V J-selected sample into two redshift subsamples, we see the possible redshift evolution of its amplitude at z ≥ 3.
The effective radii at M ⋆ = 5×10 10 M ⊙ from the bestfit are ∼ 0.6 kpc, ∼ 0.2 dex smaller than that of z ∼ 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 U V J-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).Firstly, 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) shows that the median size of green valley galaxies lies between star-forming and quiescent galaxies.If this trend persists at z ≥ 3, such larger galaxies flatten the size-mass relation.Secondly, this can be due to the possible redshift dependence of the stellar mass distribution of the sample, i.e., less low-mass galaxies  2021) in 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 ≥ 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).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.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 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 van der Wel et al. ( 2014) since the slope is still flatter when the sample is split into two redshift subsamples (Figure 7), this increase the size at the low stellar mass compared The small colored circles are the individual quiescent galaxies selected by U V J or N U V U V J.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⋆/M⊙) > 10.3, log (M⋆/M⊙) > 11.25, 10.8 < log (M⋆/M⊙) < 11.3, and 11.0 < log (M⋆/M⊙) < 11.31 for van der Wel et al. (2012), Marsan et al. (2019), Lustig et al. (2021), andEsdaile 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).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 U V J and N U V U V J-selected quiescent galaxies, respectively.
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 so as Nedkova et al. (2021), whereas van der Wel et al. (2014) apply the Axis ratio in F277W filter as 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 U V J-selected and N U V U V Jselected 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;Lustig et al. 2021;Esdaile et al. 2021).
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.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 γ.
The high median Sérsic index (n ∼ 3 − 4) and the axis ratios (q ∼ 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 U V J-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.65  −0.91 for SXDS-27434 (Appendix D) and n = 2.3 ± 0.3 for ones 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 JW ST enables us to explore the structural properties of quiescent galaxies z > 3.This study serves as a first look at this topic with JW ST .We only focus on the half-light radii here, but note that future work based on the JW ST 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 ≥ 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 ≥ 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.
We appreciate the anonymous referee for helpful comments and suggestions that improved the manuscript.(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(Perrin et al. , 2014) )  Stars used to compare with the model are selected based on the following four criteria.(1).Objects with the signalto-noise ratio of 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 FLUX RADIUS from sep.We define the range of the sequence as the 1σ 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 half-width at half maximum 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 signal-to-noise cut for the magnitude.At that wavelength range, the PSFs are determined homogeneously among fields.
B. SUMMARY OF FITTING RESULT Figure 13, Figure 14, and Figure 15 show the fitting results of all 26 objects which are used in this study.F200W and F277W are selected here because they correspond to the rest-frame ∼ 0.5 µm at 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 16 shows the relative uncertainties of the effective radius, Sérsic index, and axis ratio for mock Sérsic profile 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.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 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 (R.A., Dec. = 34.29871,−4.98987) was observed with the K s filter for 2040 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 (NDIT × DIT = 6 × 10 s).We reduce the data as described in Brammer et al. (2016) 2 .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 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 FWHM = 0. ′′ 34.The best-fit model and the residuals are shown in Figure 17.We measure an effective semi-major axis of r 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× 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. ′′ 0.8 or 0. ′′ 1) and a fixed or freely varying sky level.The target is resolved, and the 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 ± 0.20 kpc).

Figure 1 .
Figure1.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 ′′ × 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.

Figure 2 .
Figure 2. Left panel: U V J diagram(Williams et al. 2009).Squares represent quiescent galaxies selected by the U V J only, whereas circles represent those selected by U V J and N U V U V J.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⊙) > 9.5 at 3 < z < 5. Right panel: N U V U V J diagram(Gould et al. 2023).The meanings of the markers are the same as in the left panel.

A
Figure 3.A Galfit fit of #7995 as an example, which is U V J-selected quiescent galaxy at z = 3.21 with log (M⋆/M⊙) = 10.54 in 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 ′′ × 4.0 ′′ .The images are scaled by the noise (±15σ) 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.

Figure 5 .
Figure5.Wavelength dependence of the effective radius (γ = ∆ log r eff /∆ log λrest) as a function of stellar mass for U V J (left) and N U V U V J (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 is the relation of quiescent galaxies at 0.8 < z < 1.0 inKawinwanichakij et al. (2021).

FAST++ 1 .
They are in good agreement with the offset of log (M ⋆,EAZY /M ⊙ ) − log (M ⋆,FAST++ /M ⊙ ) = 0.069 +0.075 −0.066 (0.067 +0.033 −0.027 ) for our U V J (N U V U V J)selected quiescent galaxies, respectively.van der Wel et al. (2014), Mowla et al. (2019b) and Nedkova et al. ( 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 (α = 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, α, and σ log R eff ) as

Figure 6 .
Figure6.Size-mass relation of U V J-selected (left) and N U V U V J-selected (right) quiescent galaxies at z ≥ 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 are 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⊙) < 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 derWel et al. (2014), respectively.

Figure 7 .
Figure 7. Size-mass relation of U V J-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.In the right panel, the best-fit of z < 3.3 subsample is shown in a dashed blue line for easy comparison.

Table 1 .a
Summary of the obtained best-fit parameters of the sizemass relation.06 ± 0.03 0.79 ± 0.07 0.14 ± 0.03 The number of galaxies used in deriving the best-fit parameters of the size-mass relation.b Their uncertainties are taken from the 16th to 84th percentile of the total chain in the MCMC fitting.cThe best-fit parameters of early-type galaxies at 2.5 < z < 3.0 reported in van derWel et al. (2014).

Figure 8 .
Figure 8. Left panel: Redshift evolution of the effective radius at M⋆ = 5 × 10 10 M⊙ in the best-fit size-mass relation.Large red and blue circles show the cases for U V J-selected and N U V U V J-selected quiescent galaxies.Small red circles are the case for redshift-sub binned U V J-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 (α).The symbols' meanings and the data points' horizontal axis values are the same as in the left panel.

Figure 9 .
Figure 9. Compilation of the size-mass relation of quiescent galaxies up to z ∼ 4. Those from this study (z ≥ 3) are shown in the solid and dashed lines for U V J and N U V U V Jselected 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.Their color represents their median redshift.Markers represent the size measurement of individual quiescent galaxies from the literature at z ≥ 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).They are colored by their redshift, too.We note those fromForrest et al. (2022) andLustig 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.

Figure 10 .
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 U V J or N U V U V J.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⋆/M⊙) > 10.3, log (M⋆/M⊙) > 11.25, 10.8 < log (M⋆/M⊙) < 11.3, and 11.0 < log (M⋆/M⊙) < 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).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 U V J and N U V U V J-selected quiescent galaxies, respectively.
Figure11.Axis ratio in F277W filter as 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 U V J-selected and N U V U V Jselected 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;Lustig et al. 2021;Esdaile et al. 2021).
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-exclusiveaccess 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".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

Figure 13 .Figure 15 .
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 ′′ × 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 √ 2.

af
± 0.035 (0.005 ± 0.005) 4.66 ± 3.52 0.232 ± 0.313 −82.0 ± 81.0 23.94 ± 0.09 0.034 ± 0The IDs are identical to those inValentino et al. (2023).bTheyare derived with eazy-py inValentino et al. (2023).cTheeffective radius at the rest-frame 0.5 µm derived in Section 3.5.dThewavelength dependence of the effective radius derived in Section 3The 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. gThough the measured effective radii are smaller than the half-width at half maximum of PSFs in F277W, they have larger effective radii than the half-width at half maximum of PSFs in the other filters at shorter wavelength.

Figure 16 .
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.
D. NEW SIZE ESTIMATE OF SXDS-27434 AT Z = 4.01

Figure 17 .
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.

Table 2 .
Summary of the median half width at half maximum of used PSFs for images of each filter among fields.
A. SELECTION OF STARS FOR SMOOTHING PSFS AND AN EXAMPLE OF SMOOTHED PSF

Table 3 .
Summary of ID a ⊙ ) b