Enhanced Subkiloparsec-scale Star Formation: Results from a JWST Size Analysis of 341 Galaxies at 5 < z < 14

We present a comprehensive search and analysis of high-redshift galaxies in a suite of nine public JWST extragalactic fields taken in Cycle 1, covering a total effective search area of ∼358arcmin2 . Through conservative (8σ) photometric selection, we identify 341 galaxies at 5 < z < 14, with 109 having spectroscopic redshift measurements from the literature, including recent JWST NIRSpec observations. Our regression analysis reveals that the rest-frame UV size–stellar mass relation follows Reff∝M*0.19±0.03 , similar to that of star-forming galaxies at z ∼ 3, but scaled down in size by ∼0.7 dex. We find a much slower rate for the average size evolution over the redshift range, R eff ∝ (1 + z)−0.4±0.2, than that derived in the literature. A fraction (∼13%) of our sample galaxies are marginally resolved even in the NIRCam imaging (≲100 pc), located at ≳1.5σ below the derived size–mass slope. These compact sources exhibit a high star formation surface density ΣSFR > 10 M ⊙ yr−1 kpc−2, a range in which only <0.01% of the local star-forming galaxy sample is found. For those with available NIRSpec data, no evidence of ongoing supermassive black hole accretion is observed. A potential explanation for the observed high [O iii]-to-Hβ ratios could be high shock velocities, likely originating within intense star-forming regions characterized by high ΣSFR. Lastly, we find that the rest-frame UV and optical sizes of our sample are comparable. Our results are consistent with these early galaxies building up their structures inside out and being yet to exhibit the strong color gradient seen at lower redshift.


Introduction
In a hierarchical Universe, dark matter starts collapsing at initial density peaks, giving rise to the underlying structure.Baryons then start accreting in the dark matter potential wells and forming stars.Depending on the initial conditions of the gas and dark matter halos, the appearance of the resulting system may differ dramatically (e.g., Mo et al. 1998;Bullock et al. 2001b).
Within this context, the size of galaxies is a fundamental and essential proxy for understanding galaxy formation and evolution.Galaxies occupy a relative narrow portion of the size and stellar mass/luminosity plane.The distribution of galaxies within this so-called size-stellar mass/luminosity relation, the average size growth rate across cosmic time, and the distribution of other structural parameters, such as the Sérsic index and axis ratio, are key diagnostics of early galaxy formation.
In the local Universe, the large statistics enabled by the Sloan Digital Sky Survey have revealed a fundamental relation of galaxy structures and sizes with mass (e.g., Kauffmann et al. 2003).For example, it has been found that local early-type and late-type galaxies follow different slopes in the size-mass plane (e.g., Shen et al. 2003;Guo et al. 2009;Simard et al. 2011;Cappellari 2013).Such differences are believed to arise from a combination of initial conditions, evolutionary paths, and environmental influences.
In contrast, observations of galaxies at higher redshifts have revealed a variety of galaxy morphologies (Conselice et al. 2004;Wuyts et al. 2011;Guo et al. 2012;Szomoru et al. 2012).These observations indicate active physical processes within and between galaxies, establishing the structural sequences seen in the local Universe.Prior to JWST, the Hubble Space Telescope (HST) pushed the frontier of the fundamental galaxy size-mass relation to z ∼ 7 (Bruce et al. 2012;Mosleh et al. 2012;van der Wel et al. 2012van der Wel et al. , 2014;;Morishita et al. 2014;Allen et al. 2017;Yang et al. 2021).Beyond that redshift, however, the investigation has been severely limited by spatial resolution (∼1 kpc), as well as the limited number of infrared filters at >2 μm, which are critical for robustly inferring galaxy stellar masses.As such, the effort beyond that redshift has been largely limited to small samples of relatively luminous galaxies (Oesch et al. 2010;Ono et al. 2013;Holwerda et al. 2015) or a small volume through a strong gravitational lens (Kawamata et al. 2015(Kawamata et al. , 2018;;Yang et al. 2022a).
These observable properties are believed to reflect the initial conditions of the gas and dark matter halos from which galaxies form.However, over time, such fundamental properties are susceptible to contamination through a sequence of stochastic, nonlinear physical processes, including mergers.Therefore, a detailed characterization of galaxy size becomes imperative, as this may offer insights into not only the physical mechanisms at work, but also their interplay with the interstellar medium (ISM; e.g., Marshall et al. 2022;Roper et al. 2022) and the nature of dark matter (e.g., Shen et al. 2024).
Early results from Cycle 1 have already demonstrated the remarkable capabilities of JWST, with its red sensitivity and resolution, revealing early galaxy morphologies down to scales of 100 pc, throughout rest-frame UV to optical wavelengths (e.g., Finkelstein et al. 2022;Naidu et al. 2022;Yang et al. 2022b;Huertas-Company et al. 2023;Morishita & Stiavelli 2023;Robertson et al. 2023;Tacchella et al. 2023;Treu et al. 2023).Given the substantial number of observations completed in Cycle 1, this study aims to undertake a comprehensive analysis of galaxy sizes, offering the first large-scale systematic study of the galaxy size-mass relation at z > 5. To achieve this, we perform an analysis based on consistently reduced data from several publicly available extragalactic fields observed during Cycle 1, encompassing a total effective area of ∼358 arcmin 2 .This extensive coverage enables us to construct a robust sample comprising 341 galaxies within the redshift range of 5 < z < 14.
The paper is structured as follows: we present our data reduction in Section 2, followed by our photometric analyses in Section 3. We then characterize the structure of the identified galaxies and infer the distributions of their structural parameters in Section 4. We investigate the inferred physical properties and discuss the origin of early galaxies in comparison with lower-z galaxies in Section 5. We summarize our key conclusions in Section 6.Where relevant, we adopt the AB magnitude system (Oke & Gunn 1983;Fukugita et al. 1996), cosmological parameters of Ω m = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 , and the Chabrier (2003) initial mass function.Distances are in proper units, unless otherwise stated.

Data
We base our analysis on nine public deep fields from JWST Cycle 1.For all fields, except for the GLASS/Ultradeep NIRSpec and NIRCam ObserVations before the Epoch of Reionization (UNCOVER) and JWST Advanced Deep Extragalactic Survey (JADES)-GOODS (GDS) fields, where the final mosaic images are made publicly available by the teams, we retrieve the raw-level images from the MAST archive and reduce those with the official JWST pipeline, with several customized steps as detailed below.We then apply our photometric pipeline, borgpipe (Morishita 2021), on all mosaic images to consistently extract sources and measure fluxes.Our final high-z source candidates are selected by applying the Lyman-break dropout technique (Steidel et al. 1998), supplemented by photometric redshift selection, as implemented by Morishita & Stiavelli (2023).

Uniform Data Reduction of NIRCam Images
In each field, we start with raw (uncal.fits)images retrieved from MAST.We use Grizli (ver 1.8.3) to reduce the raw images to generate calibrated images (cal.fits).In this step, in addition to the official pipeline's DETECTOR step, Grizli includes additional processes for flagging artifacts, such as snowball halos, claws, and wisps.
We then apply bbpn14 on the calibrated images, to subtract 1/f noise.The tool follows the procedure proposed by Schlawin et al. (2021).Briefly, bbpn first creates object masks; then it calculates the background level in each of the four detector segments (each corresponds to the detector amplifiers) and subtracts the estimated background; it then runs through the detector in the vertical direction and again subtracts the background estimated in each row (which consists of 2048 pixels minus masked pixels); last, to compensate for any local oversubtraction of sky near, e.g., bright stars or large foreground galaxies, bbpn estimates a spatially varying background and subtracts it from the entire image.
After the 1/f-noise-subtraction step, we align the calibrated images using the tweakreg function of the JWST pipeline.For large-mosaic fields (i.e., those of Public Release Imaging for Extragalactic Research, or PRIMER, and the Cosmic Evolution Early Release Science, or CEERS, Survey), we divide the images into subgroups beforehand and process each of those separately, to optimize computing speed and memory usage.In those fields, images are split into subgroups based on the distance of each image to the other images.We here set a maximum distance of 6¢ for images to be in the same subgroup.We ensure that the images taken in the same visit (i.e., eight detector images for the blue channel and two detector images for the red channel of NIRCam) are grouped together, as their relative distance should remain consistent in the following alignment step.
The images in each subgroup are aligned on a filter-by-filter basis.We provide tweakreg with a set of images associated with the source catalogs, generated by running SExtractor (Bertin & Arnouts 1996) on each image.This enables us to eliminate potential artifacts (such as stellar spikes and saturated stars), which are often included by the automated algorithm in tweakreg, and to secure the alignment calculation by using only reliable sources.For all subgroups, each image is aligned to, when available, the WCS reference of a single, contiguous source catalog taken from large-field-of-view ground-based imaging (see below).This is to avoid alignment issues in some fields and/or subregions, e.g., misalignment in overlapping regions caused by insufficient reference stars.It is noted that tweakreg estimates a single alignment solution for images that are taken in the same visit and applies it coherently to those images, such that the distance between the imaging detectors remains the same for all visits.
Once the images are aligned to the global WCS reference, we drizzle and combine the images into a (sub)mosaic using the pipeline step IMAGE3.The pixel scale and pixel fraction for drizzling are set to 0.0315 and 0.7 for all filters.For the fields that have multiple subgroups, we then create a single mosaic using the Python function reproject.Last, to eliminate any residual shifts, we once again apply tweakreg to a set of multiband mosaics, but based on the source catalog generated using the F444W mosaic.The images are resampled after the final alignment in the same pixel grid using reproject.

JWST NIRCam Extragalactic Fields
To ensure our selection of high-redshift sources is as consistent as possible, we consider fields that have images in at least six filters (F115W, F150W, F200W, F277W, F356W, and F444W).Some fields have additional blue filters (F070W and F090W) and several medium bands (F300M, F335M, F410M, F430M, and F480M), which extend the search range toward lower redshift and improve photometric redshift estimates.When spectroscopic redshift measurements are available (from either ground or recent JWST observations), we include and use them for photometric flux calibration (Section 3.1), as well as for sample selection (spec-z supersede dropout or photo-z).

PAR1199
The PAR1199 field (11:49:47.31,+22:29:32.1)was taken as part of a Cycle 1 Guaranteed Time Observations (GTO) program (PID 1199; Stiavelli et al. 2023) in 2023 May and June, attached as coordinated parallel to the NIRSpec primary, and released immediately without a proprietary period.The NIRCam imaging consists of eight filters, including F090W and F410M, with a total science exposure of ∼16 hr with the MEDIUM8 readout mode.Due to scheduling constraints, the parallel field falls in a field where only shallow (1-2 orbits) HST Advanced Camera for Surveys F606W and F814W images are available.In addition, because of a bright star that is located near the edge of Module A, the sensitivity limit of the blue channel is slightly shallower in this module compared to the other Module, due to the increased scatter light and artifacts.The impact is found to be less in red filters.Galactic reddening is E(B − V ) = 0.024 mag (Schlegel et al. 1998).The images are aligned to bright sources in the Pan-STARRS Data Release 2 (DR2) catalog (Chambers et al. 2016).The effective area in the detection image is 10.1 arcmin 2 .

J1235
The J1235 field (12:35:54.4631,+04:56:8.50) is one of the low-ecliptic-latitude fields that were observed as part of a commissioning program used for NIRCam flat field (PID 1063; PI: Sunnquist).Among the fields available from this program, the J1235 field offers deep NIRCam coverage by 10 filters, including F070W, F090W, F300M, and F480M.Two separate visits were made in 2022 March-April and May.During the first visit, the telescope mirror alignment was not complete, and thus we exclude the data taken during that visit.With the second visit, the exposure time goes as deep as 5.8 hr in a single filter (with 50.9 hr in total for the entire field), making it one of the deepest NIRCam multiband fields.Galactic reddening is E(B − V ) = 0.024 mag.
The total field coverage extends to 34 arcmin 2 ~in the effective detection area, about ∼3.7 NIRCam footprints; however, some filter images are shifted from others, making the effective area for dropout selection dependent on the target redshift.The NIRCam images were aligned to bright sources in the Pan-STARRS DR2 catalog.

North Ecliptic Pole Time-domain Field
The NIRCam imaging in the North Ecliptic Pole Timedomain field (17:22:47.896, +65:49:21.54;Jansen & Windhorst 2018) was taken as part of a GTO program (PID 2738; Windhorst et al. 2023).The imaging data used here were taken and immediately released after the first epoch of the visit, consisting of eight filters, including F090W and F410M.For the WCS alignment of our JWST data, we use a publicly available catalog that consists of sources observed with the Subaru/Hyper Suprime-Cam instrument (Miyazaki et al. 2018) as part of the Hawaii EROsita Ecliptic pole Survey (PIs: G. Hasinger & E. Hu). 15 Galactic reddening is E(B − V ) ∼ 0.028 mag.The effective area in the detection image is 16.9 arcmin 2 .

Primer-UDS and COSMOS
Two large NIRCam mosaics are scheduled in a Cycle 1 General Observer program, PRIMER (PID 1837; PI: Dunlop).PRIMER observed two extragalactic fields, the CANDELS UKIDSS Ultra-Deep Survey (UDS) and COSMOS fields (Grogin et al. 2011;Koekemoer et al. 2011).The visits are configured with a consistent set of filters (eight filters, including F090W and F410M) and exposure time.In this study, we use the data taken during the first visit in both fields, which cover most of the entire planned fields.However, several images were identified as failed guide star acquisitions and were removed from our reduction.
The UDS mosaics are aligned to the SPLAXH SXDS catalog (Mehta et al. 2018).We include the spec-z measurements available in the same catalog.The COSMOS mosaics are aligned to the COSMOS 2020 catalog (Weaver et al. 2023).We include spec-z measurements made by the Keck/DEIMOS instrument and published in Hasinger et al. (2018).Galactic reddening is E(B − V ) ∼ 0.023 and 0.017 mag, respectively.The effective areas in the detection images are 128.2arcmin 2 and 93.9 arcmin 2 .

Next Generation Deep Field
The Next Generation Deep Extragalactic Exploratory Public (NGDEEP) survey (PID 2079; PI: Finkelstein; Bagley et al. 2023b) is a deep spectroscopic + imaging program using the NIRISS WFSS as the primary mode and the NIRCam imaging attached as coordinated parallel in the HUDF-Par2 field (Stiavelli 2005;Oesch et al. 2007;Illingworth 2009).In this study, we use the epoch 1 NIRCam imaging data, whereas the epoch 2 imaging is currently scheduled for early 2024.
The NIRCam field consists of six filters.Due to the use of the DEEP8 readout mode for many of the NIRCam exposures, a small portion of the final images are severely contaminated by cosmic-ray hits, which moderately reduces the effective field area and increases the contamination fraction in the high-z source selection.We include spec-z measurements made by the VANDELS collaboration (Pentericci et al. 2018).Galactic reddening is E(B − V ) ∼ 0.007 mag.The effective area in the detection image is 10.2 arcmin 2 .

CEERS
The CEERS Survey (Bagley et al. 2023a;Finkelstein et al. 2023) is an Early Release Science program (PID 1345; PI: Finkelstein).The data set consists of eight NIRCam filters in the Extended Groth Strip field, previously studied with HST, including as part of CANDELS.The images are aligned to the WCS of the HST F606W image released by the CEERS team (HDR1), which was originally aligned to the GAIA Early Data Release 3 WCS.
We include spec-z measurements from multiple studies (Skelton et al. 2014;Momcheva et al. 2016;Roberts-Borsani et al. 2016;Larson et al. 2022) as well as recent JWST spectroscopy studies (Arrabal Haro et al. 2023a, 2023b;Fujimoto et al. 2023;Harikane et al. 2024Harikane et al. , 2023;;Kocevski et al. 2023;Larson et al. 2023;Nakajima et al. 2023;Sanders et al. 2023;Tang et al. 2023).The effective area in the detection image is 117.5 arcmin 2 .Borsani et al. 2022c).In this study, we utilize the public imaging data made available by the GLASS-JWST team (Merlin et al. 2022;Paris et al. 2023).Their reduction processes include several customized steps, to eliminate detector artifacts.The data set consists of eight NIRCam filters, including F090W and F410M.We use the public lens model by Bergamini et al. (2023a) to correct lens magnification of the background sources.We include spectroscopic measurements made available in the literature (Braglia et al. 2009;Owers et al. 2011;Schmidt et al. 2014;Richard et al. 2021;Bergamini et al. 2023b), including those from recent JWST observations (Roberts-Borsani et al. 2022b, 2022c;Jones et al. 2023;Mascia et al. 2023;Morishita et al. 2023).Galactic reddening is E(B − V ) ∼ 0.013 mag.The effective area in the detection image is 48.4 arcmin 2 .We hereafter refer to the field as A2744, for the sake of simplicity.

JADES-GDS
We include a deep field from JADES (Eisenstein et al. 2023;Robertson et al. 2023;Tacchella et al. 2023).As of the time of writing, NIRCam imaging data in one of the deep fields in the GOODS-South field (3:32:39.3,−27:46:59) are publicly available (Hainline et al. 2023;Rieke 2023).We retrieve the fully processed images and spectroscopic catalogs made available by the team.The data set consists of nine NIRCam filters, including F090W, F335M, and F410M.We include spectroscopic sources listed in the microshutter assembly (MSA) spectroscopic catalog (Bunker et al. 2023), as well as those from the VANDELS survey (Pentericci et al. 2018).Galactic reddening is E(B − V ) ∼ 0.007 mag.The effective area in the detection image is 26.7 arcmin 2 .

Photometry
The photometric catalog in each field is constructed following Morishita & Stiavelli (2023), using borgpipe (Morishita 2021).Briefly, we first prepare a detection image for each field by stacking the F277W, F356W, and F444W filters weighted by each of their rms maps.Source identification is made in the detection image using SExtractor (Bertin & Arnouts 1996).Fluxes are then estimated for the detected sources with an r = 0 16 aperture.For the aperture flux measurement, images are point-spread function (PSF)-matched to the PSF size of F444W beforehand.The image convolution kernels are generated by using pypher (Boucaud et al. 2016) on the PSF images generated by webbpsf (Perrin et al. 2014).
We follow the standard procedure used in the literature (Trenti et al. 2012a;Bradley et al. 2012;Calvi et al. 2016;Morishita et al. 2018;Morishita & Stiavelli 2023), including correction for Galactic extinction and rms scaling that accounts for correlation noise in drizzled images.The limiting magnitudes of the images are measured by using the same aperture size and reported in Table 1.The aperture fluxes of individual sources are then corrected by applying the correction factor C = f auto,F444W /f aper,F444W universally to all filters, where f auto,F444W is the FLUX_AUTO of SExtractor, measured for individual sources.With this approach, the colors remain those measured in apertures, whereas the total measurements derived in the following analyses (such as stellar mass and star formation rate) represent the total fluxes (see also Section 3.3).Last, since several fields have a number of spectroscopic objects across a wide redshift range, we run eazypy, a Python wrapper of the photometric redshift code eazy (Brammer et al. 2008), to fine-tune the fluxes across all filters.A set of correction factors for all filters is derived in each field, from the redshift fitting results for those with spectroscopic redshift.The derived correction factors are found to be <2% relative to the pivot filter, here set to F150W in all fields, only requiring minor correction.

Selection of High-redshift Galaxy Candidates
In what follows, we present our selection of high-redshift galaxies and galaxy candidates.To construct a robust photometric sample, we adopt a twofold selection method, which has been established in our previous studies (Morishita et al. 2018;Ishikawa et al. 2022;Roberts-Borsani et al. 2022a) and is described in the following two subsections.

Lyman-break Dropout Selection
Here we identify dropout sources in four redshift ranges.For those detected at signal-to-noise ratio (S/N) >4 in the detection band, we apply one of the following criteria.
F070W-dropouts (5.0 < z < 7.2): 2 .0 , 090, 070 < z 6; set = F115W-dropouts (9.7 < z < 13.0): 2 .0 , 150, 115, 090, 070 < z 10; set = where the S/N is measured in an r = 0 16 aperture.In each redshift range, we ensure secure selection by requiring an S/N > 8 detection in a filter that covers the rest-frame UV (∼1600 Å, but not including the blue side of the Lyman break).This stringent requirement ensures high completeness (>50%) and reliable size measurements (Appendix C).For the source to be selected as dropout, we require the nondetection of fluxes (S/N < 2) in all available dropout filters (listed as subscripts above).Furthermore, to secure the nondetection, we repeat the nondetection step with a smaller aperture, r = 0 08 (∼2.5 pixels).Note that a photometric selection is not attempted for fields where no dropout filter is available (but see Section 3.2.3).
A major difference from the conventional Lyman-break technique in the literature (e.g., Bouwens et al. 2023) is that our selection method above does not cut samples based on the color of the rest-frame UV, but only on the strength of the Lyman break.The choice is made to preserve as many potential sources as possible and make the selection comprehensive-for example, in a conventional color-cut selection, sources may be dismissed for their color that barely miss the selection window, even if the color is consistent within the photometric uncertainty.This also means that the fraction of low-z interlopers misidentified as high-z sources is likely increased with regard to the standard technique.Therefore, we further secure the sample in the following step.

Photometric Redshift Selection
We here secure the dropout sources by applying photometric redshift selection to the dropout sources selected above.This is to minimize the fraction of low-z interlopers, such as galaxies of relatively old stellar populations (e.g., Oesch et al. 2016) and with dust extinction or foreground dwarfs (e.g., Morishita et al. 2020).Such interlopers are often distinguished by their distinctive red color, readily discernible in our wavelength coverage with NIRCam.
To estimate photometric redshifts, we run eazy with the default magnitude prior (Figure 4 in Brammer et al. 2008).The fitting redshift range is set to 0 < z < 30, with a step size of ( ) z log 1 0.01 + = . By comparing photometric redshifts with spectroscopic ones, we find that the template library provided by Hainline et al. (2023) offers an improved photometric redshift accuracy over the default (v1.3) template library, and thus in this work we adopt the former.
Following Morishita et al. (2018), we exclude sources that satisfy p(z < z set ) > 0.2, i.e., the total redshift probability at z < z set is greater than 20%, where z set is set separately for each redshift selection range (see Section 3.2.1).To eliminate potential contamination by cool (T-/L-/M-type) stars (i.e., brown dwarfs), we follow Morishita (2021) and repeat the phot-z analysis with dwarf templates.A set of dwarf templates taken from the IRTF spectral library (Rayner et al. 2003) is provided to eazy and fit to the data with redshift fixed to 0. The fitting result is inferred for every photometric source that is unresolved (see Section 4.4), and the source is removed if the χ 2 /ν value is smaller than the one from the galaxy template fitting above.We have excluded 60 sources in this step.
Last, we visually inspect all sources that pass the two selections above.In this step, we exclude any suspicious sources whose flux measurements may be significantly affected, including those with residuals of CRs, close to/ overlapping with a brighter galaxy (caused by deblending), misidentified stellar spikes, and those near the detector edge where a part of the source is truncated.We have discarded 342 sources.

Spectroscopic Sample
In addition to the photometric sample constructed above, we include those with spectroscopic redshift confirmed by previous spectroscopic observations, as described in Section 2.2.We add sources when their spectroscopic redshift is within the redshift range defined for each selection window and when they satisfy detection criteria in the detection (S/N >4) and rest-frame UV (>8) bands.
The addition of spec-z sources aids in particular the F070Wdropout sample, which would need F070W as a nondetection filter.All fields except for J1235 do not have the filter coverage, leaving the sample size relatively small without spectroscopically confirmed objects.On the other hand, this could introduce a potential bias toward strong-line emitters.However, in Section 4.2, we investigate this in our size-mass analysis and find that the addition of spec-z sources does not impact any of our final conclusions.

Size Measurements
Our primary analysis is based on the size measurement of galaxies.Following standard practice, we adopt the Sérsic profile (Sérsic 1963): where the size is characterized by the effective radius, R e , which encloses half of the total light of the galaxy, n is the Sérsic index, and b n is an n-dependent normalization parameter.We model the 2D light profile of each galaxy using galfit (Peng et al. 2002(Peng et al. , 2010)).We follow Morishita et al. (2014Morishita et al. ( , 2017) ) for detailed procedures, with a few modifications to accommodate efficiency and accuracy.Briefly, for each galaxy, we first generate image cutouts (here set to 151 × 151 pixels in size, equivalent to 4 8 × 4 8) of the original (i.e., pre-PSFmatched) science map, rms map, and segmentation map.We fix the Sérsic index n to 1, a value that is found to offer a reasonable fit to high-z Lyman-break galaxies in the literature (Shibuya et al. 2015;Ono et al. 2023a;Yang et al. 2022a).As a test, we repeated the analysis with n as a free parameter and indeed found its distribution to be centered around ∼1.However, this led to an increased fraction (∼13%) of unsuccessful fits, where the solutions either did not converge or converged to unrealistic parameters (i.e., n < 0.2 or n > 10).To mitigate potential bias from this constraint, we reevaluate the uncertainty of each size measurement by adding the difference in R e resulting from the two procedures (n fixed and n free) in quadrature.Consequently, those with n deviating significantly from 1 incur a larger uncertainty in the size measurement.
The PSF image generated by webbpsf of the corresponding filter is fed to galfit, for convolution of the model profile at each iteration.The PSF image is generated for each field by retrieving the Optical Path Difference files of the observed date. 16As reported in several studies (Ono et al. 2023a;Ito et al. 2023;Tacchella et al. 2023), we find that the default output of webbpsf exhibits a narrower PSF profile compared to observed stars in our reduced images.This is potentially caused by, e.g., drizzling/resampling of the actual images, as well as jitter in pointing, which could have a more significant effect in a longer exposure.We therefore find an optimal PSF that describes the actual PSF size of our image by tweaking the jitter_sigma parameter of webbpsf.To do this, we visually identify unsaturated, bright stars that do not have any companion within <50 pixels.We then fit these stars with various PSF models by using galfit and select a model that offers the minimum χ 2 value.While it is ideal to repeat the analysis and determine the jitter value in each field, some fields do not have sufficient numbers of stars that can be used for this.We thus adopt the median value that is determined by the stars in all fields for each filter.Figure 1 shows an example radial flux profile of an actual star, compared with those of webbpsf generated with different σ jitter values.The final jitter_sigma value is set to ∼0 022 for the blue channel filters and ∼0 034 for the redder channel.
Neighboring sources that are close to and relatively bright compared to the main galaxy are fit simultaneously, while the rest of the sources in the stamp are masked using the segmentation map generated by SExtractor above.We include any neighboring source at distance d from the main galaxy when its flux is above the limiting flux, defined as nei,lim main s nei g = a with γ = 0.8, α nei = 2.0, and d s = 60 pixels.
We run galfit in the order of the target source magnitudes.The fitting results of the primary galaxy are continuously stored, so that the parameters for the repeated galaxies are fixed to the previously determined values when they appear in the cutout of a fainter galaxy later in the fitting session.
For each galaxy, we repeat the fit in two filters that correspond to the rest-frame UV and optical wavelengths.We then inspect all fitting results to ensure the measurements.We have flagged 24 sources that show significant residuals, e.g., from multiple clumps within the defined segment region and/ or clear features of interaction with nearby sources.These flagged sources are excluded from statistical analyses in the following sections.The measured sizes are presented in Appendix B.

Physical Properties Inferred by Spectral Energy Distribution Analysis
We infer the spectral energy distribution (SED) of the individual galaxies through SED fitting using photometric data that cover 0.6-5 μm.We use the SED fitting code gsf (ver1.8;Morishita et al. 2019), which allows flexible determinations of the SED by adopting binned star formation histories, also known as nonparametric.gsf determines an optimal combination of stellar and ISM templates among the template library.For this study, we generate templates of different ages, [10,30,100,300,1000,3000] Myr, and metallicities at an increment of 0.1 by using fsps (Conroy et al. 2009;Foreman-Mackey et al. 2014).A nebular component (emission lines and continuum) that is characterized by an ionization parameter [ ] U log 3, 1 Î -is also generated by fsps (see also Byler et al. 2017) and added to the template after multiplication by an amplitude parameter.The dust attenuation and metallicity of the stellar templates are treated as free parameters during the fit, whereas the metallicity of the nebular component is synchronized with the metallicity of the stellar component during the fitting process.
The posterior distribution of the parameters is sampled by using emcee for 10 4 iterations, with the number of walkers set to 100.The final posterior is collected after excluding the first half of the realizations (known as burn-in).The resulting physical parameters (such as stellar mass, star formation rate, rest-frame UV slope β UV , metallicity, dust attenuation A V , and mass-weighted age) are quoted as the median of the posterior distribution, along with uncertainties measured from the 16th to 84th percentile range.The star formation rate of individual galaxies is calculated with the rest-frame UV luminosity (∼1600 Å) using the posterior SED.The UV luminosity is corrected for dust attenuation using the β UV slope, which is measured by using the posterior SED, as in Smit et al. (2016): The attenuation-corrected UV luminosity is then converted to the star formation rate via the relation in Kennicutt (1998): Last, we correct both stellar mass and star formation rate measurements to the total model magnitude derived by galfit, as in Morishita et al. (2014), by multiplying the correction factor:

Overview of the Final Sample
In Table 2, we report the number of final sources selected in each field and dropout selection window.From all fields, we identify 101 F070W-dropout sources (81 of which are spec-zconfirmed), 196 F090W-dropout sources (22), 42 F115Wdropout sources (5), and 2 F150W-dropout sources (1).As we present in the following sections, the sample spans a wide range of stellar mass ( * M M 6.8 log 10.4  


) and absolute UV magnitude (−23  M UV / mag  − 17).The final sources are identical when a larger aperture (0 32) is adopted for the source selection in Section 3.2.1.
In Figure 2, we show the distribution of the final sample in the stellar mass-star formation rate plane.Despite the wide mass (more than 3 orders of magnitude) and redshift (5 < z < 14) ranges, our galaxies are found to distribute along a sequence, suggesting that most of our sample galaxies are a typical star-forming population.To compare the locations of our galaxies with those at lower redshift, we derive a linear regression, with the slope fixed to 0.81, i.e., the one for z = 5

Table 2
Numbers of the Final Sources in Fields
While the full details of the selected sources will be presented in a forthcoming paper, we highlight two F150W-dropout galaxies (z  13).One, JADESGDS-30934, is spectroscopically confirmed to be at z = 13.2 (Curtis-Lake et al. 2023).The other object, PRIMERCOS-38203, is a newly identified photometric candidate source at z 13.8 1.0 in the PRIMER-COSMOS field (Figure 3).Despite the relatively shallow depth in the field, the source exhibits clean nondetection in F090W, F115W, and F150W and high-S/N detection in F200W (S/N = 9.2) and in the IR detection band (S/N = 16.1).The observed UV magnitude (M UV = − 20 mag) and the derived stellar mass ) are both moderate and comparable to other galaxy candidates at these redshifts (e.g., Bouwens et al. 2023;Finkelstein et al. 2023;Morishita & Stiavelli 2023).While the F200W-dropout selection covers up to z ∼ 18.6, none is identified at z > 14 in our selection.The number density estimates of the identified sources will be presented in a forthcoming paper.Another potentially interesting source is PRIMERUDS-121885 at z 10.9 0. .This object has relatively red F356W −F444W color (0.76 mag), which implies the presence of old populations.However, we caution that the source is identified in PRIMERUDS, which is relatively shallow among the fields.
In addition, at the redshift of the source, the F444W flux could also be attributed to strong Hβ+[O III] emissions, which would lead to a smaller stellar mass.
Last, we observe a mild enhancement of sources at z ≈ 7-7.6.This is partially attributed to strong Hβ+[O III] emitters being more sensitive to our selection, due to the two medium-band filters (F410M and F430M), by their making the photometric redshift relatively more constrained.Besides, there is an overdensity of emitters identified in the same redshift range in one of the fields (K.Daikuhara 2023, in preparation.).

Size-Stellar Mass Distribution of Galaxies at 5 < z < 14
In Figure 4, we show the distribution of galaxies in the sizemass plane for the four redshift ranges.In each redshift panel, we show the size measured in the filter that corresponds to restframe ∼1600 Å i.e., F115W for the F070W-dropout, F150W for the F090W-dropout, F200W for the F115W-dropout, and F277W for the F200W-dropout selection.We adopt the effective radius measured along the major axis, to mitigate the effect by inclination.
The measured size spans a broad range, R log kpc 2 e ~to 0.3.Remarkably, at * M M log 9  < , many galaxies are characterized by R e < 0.3 kpc (<0 07), which is below the resolution limit afforded by HST/WFC3-IR.Thus, the spatial resolution of NIRCam, ∼0.11-0.14kpc, is essential to study typical star-forming galaxies at these redshifts (see also Section 4.4).
We investigate their distribution on the size-mass plane by linear regression analyses.By following Shen et al. (2003;also Ferguson et al. 2004;van der Wel et al. 2014), we parameterize size by a log-normal distribution as a function of stellar mass and redshift: where we describe the intercept as We adopt a pivot mass M 0 = 10 8 M e .The model distribution, ( ( ) , prescribes the probability distribution for observing R log e for a galaxy with the stellar mass M * with an intrinsic scatter of logR e s .By making the intercept B(z) a function of redshift, we are able to evaluate the size distribution with a consistent slope determined by the entire sample.Although the redshift evolution of the slope is of interest, previous studies have reported that little changes over a much broader cosmic time range than what is explored in this study (e.g., van der Wel et al. 2014;Shibuya et al. 2015).As we demonstrate below, our findings also reveal no significant evolution in slope from these studies, supporting our employing a single slope across the entire redshift range.
The regression is determined by using emcee (Foreman-Mackey et al. 2014), with the number of walkers being set to N = 50 and the iterations to 10 4 .We only include those not flagged in the structural fitting analysis in Section 3.3.Measurement uncertainties in size, stellar mass, and redshift are used in the calculation of likelihood.The derived regression is shown in Figure 4, along with the measured sizes in four redshift panels.In Figure 5, we also show the redshift-corrected size, ( ) R z 1 e z + a , as this represents the actual variants evaluated in the fitting.We report the determined parameters in Table 3.
The derived slope, ∼0.20 ± 0.03, is similar to what was found in Mosleh et al. (2011;∼0.3)for Lyman-break galaxies at z ∼ 3.5 and van der Wel et al. (2014; 0.18-0.25)for late-type galaxies at z ∼ 0.3-3, despite the latter being observed in different rest-frame wavelengths (but see the following and Section 5.3).The slope of the size-mass relation would reflect the interplay between the intrinsic compactness and concentrated dust attenuation in the core of massive star-forming galaxies at high redshift (Roper et al. 2022).In fact, negative slopes of the size-mass and size-luminosity relations have been found in cosmological simulations, e.g., the BlueTides (Marshall et al. 2022), IllustrisTNG (Popping et al. 2022;Costantin et al. 2023), FLARES (Roper et al. 2022), and THESAN (X.Shen et al. 2023, in preparation) simulations.Our finding of a positive slope is consistent with the idea that some massive galaxies in our sample may already possess a moderate amount of dust in their cores.We note that the derived intrinsic scatter, , is relatively large (∼0.3) compared to those at lower redshift in the literature (∼0.2).The scatter in size distribution ought to reflect the initial conditions of the dark matter halos, such as the distribution of the spin parameter (e.g., Bullock et al. 2001a).Thus, the observed larger scatter would imply the presence of galaxies that experienced nonlinear behavior, such as mergers and other dissipative processes, and deviate from the distribution predicted by a simple galactic disk formation model (e.g., Mo et al. 1998).In fact, we observed a number of galaxies that are barely resolved in the NIRCam images, some of which are located far below (>1.5σ) the derived regression (Section 4.4).We repeat the regression analysis by excluding these unresolved sources, but with the slope fixed to the one derived above.The derived scatter from this regression is ∼0.25, moderately reduced from the analysis of the full sample (Table 3).
Our sample includes both spectroscopic and photometric sources.In particular, the majority (∼90%) of the F070W-d sample is spectroscopic, due to the lack of F070W filter coverage.To investigate the impact of the spectroscopic sources, we repeat the regression analysis by only using photometric sources.We find a consistent result (α = 0.23 ± 0.04, β z = −0.10 ± 0.32, and α z = −0.56 ± 0.33) within the uncertainty range.The slight decrease in α z (+0.14; stronger z-evolution) can be attributed to the fact that the spectroscopic sources are smaller in size than the photometric sources and they dominate the lowest-redshift bin.These differences in the regression result do not change our conclusion.

Redshift Evolution of Galaxy Sizes
In Figure 6, we show the redshift trend of the UV sizes derived through the regression analysis above.We note that the redshift evolution represents, by design, for the mass-corrected size, a .We find that the derived redshift evolution, ( ) z 1 z µ + a with α z ∼ −0.4 ± 0.2, is much less significant than the one derived for the rest-frame UV (∼2100 Å) size of Lyman-break galaxies at 0 < z < 7, with α z = −1.2± 0.1 (Mosleh et al. 2012; see also Oesch et al. 2010;Shibuya et al. 2015, who found a similar value).
The primary cause of the discrepancy can be attributed to the redshift range of the sample.The aforementioned studies derived the redshift evolution by including lower-redshift galaxies (z  1 in Mosleh et al. 2012;Shibuya et al. 2015; and 2 < z < 8 in Oesch et al. 2010).Indeed, Curtis-Lake et al.
(2016) found a much slower evolution, α z , for the Lymanbreak galaxy (LBG) sample at 4 < z < 8 and argued that the discrepancy is partially attributed to a stronger evolution in UV sizes at z < 5 (see also Oesch et al. 2010, who found little size evolution from z ∼ 7 to 6).
To investigate this, we derive the redshift evolution by combining the median sizes presented in Figure 6 and in Mosleh et al. (2012; for 1 < z < 7), and we do find a stronger evolution (α z ∼ −1.2).This supports the idea that the redshift evolution of the average UV size is much slower at z > 5 than lower redshifts.However, we note that the outcome is likely subject to a range of systematic factors, such as the weighting of each size measurement, the inclusion or exclusion of specific data points, and potential sample mismatches.
In Figure 6, we also show the evolution of the rest-frame optical sizes derived for low-mass ( * M M log 9  ~) late-type galaxies at 0 < z < 2 (van der Wel et al. 2014).Interestingly,   et al. 2014).We note that the size trend of both previous studies is also scaled to the same pivot mass.Median sizes (large symbols) are derived in each redshift window.Those identified as compact (Section 4.4) are shown by magenta symbols.the extrapolated sizes from the fit to our redshift range exhibit a significant offset, even after accounting for the mass difference between the two samples.Although the offset might be attributed to differences in the bands used for size measurements (i.e., rest-frame 5000 Å in van der Wel et al. 2014), we find that our galaxies, on average, do not show a significant offset between the two bands.This reinforces our earlier interpretation that size evolution could be more pronounced at z  5, primarily driven by the buildup of massive bulges and/ or outskirts, which would enhance color gradients, as seen at lower redshifts.We discuss this in more detail in Section 5.3.

Identifying Blue Compact Sources
In our size analysis above, we have identified a number of compact sources that are near the resolution limit of NIRCam.Previous studies using NIRCam data also reported a few of these compact sources (Castellano et al. 2022;Naidu et al. 2022;Ono et al. 2023a;Yang et al. 2022a).We define those that satisfy either of the following as compact: where ΔR is the difference of the measured size from the inferred size by the linear regression for the stellar mass and redshift of the source.Note that R e used here is the apparent size.
With these criteria, we find 44 compact sources from our sample (∼13%), including 15 that are spectroscopically confirmed.We then exclude five of the classified compact sources, whose SEDs are better fitted with a brown dwarf template over galaxy templates (Section 3.2).The classified compact sources are marked in Figure 4. Individual cutout images are presented in Appendix D.
In Figure 5, we show the distribution of the normalized size, . The distribution clearly shows an excess at smaller size when compared to the distribution of the remaining noncompact sources.We note that the fraction of identified compact sources is well above the number expected at 1.5σ for a normal distribution (i.e., 6.7%).The compact sources follow a similar distribution as other extended sources in physical properties, such as stellar mass, star formation rate (see also Figure 2), and β UV , except for the star formation surface density (Section 5.1).The ISM properties of the compact sources are further investigated in Section 5.2.

High Star Formation Efficiency Revealed by NIRCam Imaging
The star formation (rate) surface density, is known as an excellent proxy for inferring the current mode of star formation.In Figure 7, we show the distribution of Σ SFR of our sample.Overall, the median values of our galaxies are consistent with previous studies with HST (Oesch et al. 2010;Shibuya et al. 2015).However, we observe a large scatter along the vertical axis.In fact, we find a large fraction of galaxies of log 1.5 SFR  S , comparable to local ULIRG/starburst (Scoville et al. 2000;Dopita et al. 2002;Kennicutt & Evans 2012) and high-z submillimeter galaxies (e.g., Daddi et al. 2010).Such high-Σ SFR galaxies were not reported in previous HST studies at similar redshift (Oesch et al. 2010;Ono et al. 2013;Holwerda et al. 2015), except for a few cases in cluster lensing fields (e.g., Kawamata et al. 2015;Bouwens et al. 2022).Only ∼6% of the low-z sample is found at Σ SFR > 1 M e /yr/kpc 2 ; the majority of our sample is located above the same value.
Obviously, Σ SFR estimates are limited by imaging resolution, and thus the previous estimates for smaller galaxies identified in HST data often remained as lower limits (e.g., Morishita 2021;Fujimoto et al. 2022;Ishikawa et al. 2022).
The observations in Figure 7 also demonstrate the potential for JWST NIRCam observations to probe star formation activity at low stellar masses (M  10 8 M e ).This has so far been limited to indirect probes, such as studies based on follow-ups of gamma-ray burst afterglows to quantify the fraction of detected host galaxies (e.g., Tanvir et al. 2012;Trenti et al. 2012b;McGuire et al. 2016).Those studies hinted at the existence of the population of sources now directly probed by NIRCam.

Implications for Star Formation Efficiency
The compact sources identified in Section 4.4 dominate the upper range, log 1.5 SFR  S . We note that this does not stem from their star formation rates, which are, on average, comparable to those of more extended sources.Instead, it is attributed to their compact nature.These compact sources' physical characteristics are particularly intriguing, as negative feedback is likely more effective within confined systems.The observed high values thus imply efficient gas fueling within these compact sources, potentially facilitated by processes such as the loss of angular momentum through mergers (see also Section 5.2).
In addition, the mass dependence of Σ SFR could hint at the efficiency of star formation, as the system's mass is tightly linked to the regulation of star formation.In the right panel of Figure 7, we show the distribution of Σ SFR as a function of stellar mass.Of particular interest is the high mass range, ∼10 9 M e .In this mass range, the shocked gas remains hot (e.g., Birnboim & Dekel 2003;Dekel & Birnboim 2006;Stern et al. 2021), resulting in a reduced star formation efficiency, as seen at lower redshifts.The observed high values for our sample thus imply that our sources still hold a high efficiency within this mass regime.We calculate the median massdoubling time by t = M * /SFR and find ∼15-90 Myr for our sample.These considerably small mass-doubling times suggest that some galaxies in our sample could evolve to * M M log 11  ~by z ∼ 5, provided that the efficiency remains similar in the following ∼0.5-1Gyr.

Implications for the Growth of Galaxies
At these early cosmic times, comparing the radius-mass relationship of galaxies as a function of the age of the stellar population can reveal how galaxies grow (Figure 8).In hierarchical structure formation, galaxies grow through the merger of smaller-mass galaxies.The merger process results in the infall of stars, but also allows new star formation to occur through the shocking of the infalling gas.If galaxies grow predominantly by cold gas accretion (e.g., Dekel & Birnboim 2006), one would expect their size to evolve predominantly through secular evolution, which would be driven by dynamical friction between the stars and the accreted gas, i.e., galaxies would become smaller as they built up their stellar mass.We find two clear trends in our sample of high-z galaxies: 1.Both radius and age increase with stellar mass.Highmass galaxies harbor older (∼ several 100 Myr) stellar populations compared to lower-mass galaxies, which harbor populations of age ∼10-100 Myr. 2. The scatter in radii at fixed mass is larger by 0.2 dex for low-mass galaxies than for high-mass galaxies.
Since there is no selection effect that prevents the selection of small, high-mass galaxies, the implication of these two observational trends is that galaxies form inside out.Small, low-mass galaxies have a large scatter in their sizes, likely due to a combination of gas accretion and merger events.As they form stars, the stars get scattered due to three-body interactions, resulting in a growth in size and stellar mass.That is consistent both with the size evolution with cosmic time and the radius and stellar age evolution with mass.Although the sensitivity of the data at the present time is not adequate to constrain minormerger rates at these redshifts, future deeper surveys will help develop this hypothesis further.

Comparison with Local Lyman-break Analogs
It is illustrative to compare the luminosity surface density of these galaxies with similar objects in the local Universe (Hoopes et al. 2007;Overzier et al. 2009;Shim & Chary 2013).Although objects with such high surface densities of star formation exist at z ∼ 0, less than 0.01% of galaxies in the Sloan sample have Σ SFR > 1 M e yr −1 kpc −2 .Even at higher redshifts, the fraction remains small (∼6% at 0.3 < z < 3.5; Skelton et al. 2014).In contrast, the fraction is ∼100% in our sample.The z ∼ 0 objects have a median mass of 10 8.9 M e and span the full range of metallicity, from subsolar to supersolar.They are not particularly biased toward young ages, as determined from the strength of the Balmer break.Morphologically, we see some evidence of mergers in the high-z sample, with ∼20 % of galaxies having nearby  5, but the symbols are color-coded by the age ( t log ) derived in our SED analysis.The median age in each grid is shown by the larger symbols (squares), with the symbol size scaled by the number of galaxies in the grid.In the inset, the scatters in size measured at each mass bin are shown.The error of each scatter measurement is estimated by bootstrapping.Clearly, the radius and age of the stellar population increase, while the scatters in size decrease, with increasing stellar mass.
companions (besides the 24 galaxies flagged in the size analysis).Taken together, the implication is therefore that the high surface densities are gas-rich galaxies, with the gas being relatively concentrated toward the nucleus of the galaxy, indicating that the late stages of the merging process are driving the inflow of gas toward the nucleus of the galaxies.

Nature of the Compact Sources-NIRSpec Spectroscopic Analysis
We here investigate the nature of the compact sources identified in Section 4.2, specifically aiming to observe if they exhibit any evidence of active galactic nuclei (AGNs).Traditionally, point-source-like morphologies of UV-bright sources have been considered as evidence of AGNs.However, it has been shown that the compactness does not always indicate the presence of AGNs (e.g., Morishita et al. 2020) and vice versa (see, e.g., Matthee et al. 2023, for faint AGNs with extended features).As such, a comprehensive approach is required to answer the question.
A subset of our sample galaxies have spectroscopic coverage by NIRSpec/MSA, taken as part of CEERS, GLASS-ERS, and JADES.Following Morishita et al. (2023), we reduce the MSA spectra using msaexp.17For the extracted 1D spectrum of each source, we fit the line profiles of Hβ and [O III] doublets with a Gaussian, after subtracting the underlying continuum spectrum inferred by gsf in Section 3.4.The total flux of each line is estimated by integrating the flux over the wavelength range of 2 × FWHM derived from the Gaussian fit.In the following analysis, we include sources with measured line S/Ns above 3 for the [O III] λ5007 line (N = 51); when Hβ is not detected above the same significance, we quote the 3σ flux limit measured at the wavelength of the line over the same line width derived for [O III] λ5007 .
In Figure 9, we show the sample in the mass-excitation (MEx) diagram, a conventional diagnostic for AGNs and starforming galaxies (Juneau et al. 2011).Among the 51 objects, we find eight sources located within the MEx AGN region defined by Juneau et al. (2014).One of the sources, CEERS7-18822, was previously reported to have a tentative (∼2.5σ) broad component in Hβ (Larson et al. 2023), agreeing with the classification here.CEERS7-18822 exhibits extended structures in F115W, which confirms our classification of this source not to be compact.
None of our compact sources are classified as MEx-AGNs.While JADESGDS-18784 is located near the MEx border, with log [O III]/Hβ = 0.83 ± 0.34 (classified as a MEx-AGN within the uncertainty range), we do not confirm any features that immediately support the presence of AGNs (i.e., broadline features or high-ionization lines).None of the other compact sources show AGN signatures either, except for CEERS6-7832, which was reported to have a broad (∼2000 km s −1 ) component in its Hα line (Kocevski et al. 2023;as CEERS_1670). 18owever, this does not completely rule out the absence of AGNs.First, the discriminating power of the MEx diagram could be lower near the transition region.As shown by Juneau et al. (2014), for the range of the log [O III]/Hβ line ratios probed here (0.4), there could be 10%-30% of AGNs present even inside the MEx-star-forming region at * M M log 10  ~.For example, the aforementioned CEERS_7832 is found in the MEx-star-forming region.The accuracy of the AGN classification at a lower mass range is not known, due to the lack of data in the previous study.We also note that the observed ratio for our sample (log [O III]/Hβ 0.6) is relatively high, compared to the star-forming galaxies of a comparable mass in Juneau et al. (2014; see also Shapley et al. 2015).Such a high ratio can still be achieved by a stellar-only configuration, but requires, e.g., high electron density (Reddy et al. 2023).
We note that given the evolving ISM properties at these redshifts, there is likely a shift of the MEx boundary toward a higher line ratio, as is the case for z ∼ 0 to z ∼ 2. As such, it is still possible that any of our MEx-star-forming sources near the border may host an AGN, if not a broadline AGN.Furthermore, as has been demonstrated in the local Universe (Ho et al. 1997), it is possible to bury a low-luminosity Seyfert or LINER-type nucleus in a galaxy without detecting a component of broadline emission.Vice versa, some lowermass MEx-AGN sources here could turn to be MEx-starforming, due to the potential redshift evolution of the boundary.
Nonetheless, the absence of clear AGN evidence implies that the observed high ratios are driven by the high surface density of star formation.When we fit high values of [O III]/Hβ >3 with radiative shock models (e.g., MAPPING III; Allen et al. 2008), we find that shock velocities of a median of 500 km s −1 are required.Ratios of ∼10 can only be achieved with high shock velocities, which could be correlated with the high surface density of star formation (or strong AGN activity).Achieving the observed high star formation surface density (>100 M e /yr / kpc 2 ) is considered challenging, due to the presence of negative feedback.Such a high density is expected only in an extreme environment, where an abundance of gas is available, and/or in (post-)merging systems, where gas could rapidly fall in.By comparing with a numerical simulation, Ono et al. (2023a) found that such a compact galaxy is in a temporary compact star-forming phase triggered by recent major mergers.Roper et al. (2022) found a large fraction of blue compact (∼100-300 pc) galaxies and also found these galaxies to have little contribution from AGNs in the FLARES simulations (see also Marshall et al. 2022;X. Shen 2023, in preparation).Future work will compare the derived physical parameters from the emission lines with the measured star formation rate density.

Comparison of Sizes in Rest-frame UV and Optical
In this study, we have analyzed the sizes of galaxies at a restframe wavelength of ∼1600 Å.In Section 4.2, we found that the average sizes of our galaxies are much smaller than those predicted from the extrapolation of van der Wel et al. (2014) at the corresponding redshift.A possible explanation for the discrepancy could be the rest-frame wavelengths where the size is measured in the previous study (i.e., ∼5000 Å).We investigate this by repeating our size analysis, but in a filter that corresponds to ∼5000 Å, i.e., F356W for the F070Wdropouts and F444W for the F090W-dropouts.For the other two ranges at higher redshift, we use the reddest filter available (F444W), which corresponds to ∼4000 Å and ∼3000 Å, respectively.In Figure 10, we show the distribution of the size difference in these two wavelengths.We find that the difference of size in the two filters is negligible on average and thus conclude that the wavelength difference cannot explain the offset seen in the extrapolated size of van der Wel et al. (2014) observed in Figure 6.Instead, we speculate that the extrapolation of their relation may not persist in the redshift range far beyond their probed redshift range z < 3.
In fact, we have seen in Section 4.2 that the size evolution is much slower (α z = −0.4)than the one found in previous studies of rest-frame UV size (∼ −1.2; Mosleh et al. 2012;Shibuya et al. 2015).While a comprehensive analysis covering a wide redshift range would be necessary, we attribute the observed conflict to the buildup of complicated structures in galaxies, such as a central massive bulge and young starforming disk.At lower redshifts, we observe a more pronounced difference in sizes between different wavelengths, driven by radial color gradients within the systems (e.g., Vulcani et al. 2014).Similarly, Shibuya et al. (2015) found that at z = 1.2-2.1, the average UV size is smaller than the optical size by ∼20% in the low-mass regime (∼10 9 M e ), while the trend is reversed at the high-mass end (10 11 M e ; see also Szomoru et al. 2012;van der Wel et al. 2014).
Figure 10.Distributions of rest-frame UV-to-optical size ratio for three redshift bins.The first two panels compare the rest-frame UV to optical (∼5000 Å) sizes, whereas the last panel is at a shorter wavelength (∼3700 Å), due to the filter availability.Only resolved galaxies in the two filters are included.The F150W-dropouts are not shown, as neither of their sizes are resolved in F444W.Each distribution is fit with a Gaussian (red line), with the mean position indicated by a dashed vertical line.The difference with the extrapolated size in the rest-frame optical of van der Wel et al. (2014; observed in Figure 6) is shown in each panel (dotted vertical line); the difference is measured to be 3.2σ, 2.7σ, and 3.7σ, respectively.
On the other hand, the resemblance of galaxy UV and the optical sizes of high-redshift galaxies is not unexpected, given the rapid assembly of the stellar content (see also Yang et al. 2022b;Treu et al. 2023).For our galaxies, we have estimated the mass-doubling time to be 100 Myr in Section 5.1.This timescale is comparable to or even smaller than the star formation timescale that the UV tracer is sensitive to (Murphy et al. 2011;Flores Velázquez et al. 2021).Consequently, the stellar content detected by the UV will predominate over the total content integrated over the full star formation history, which is probed by observations at optical wavelengths.The majority of our galaxies are actively building their structures inside out, developing rapidly enough to maintain coherence.
Last, we investigate the rest-frame optical size of the compact sources identified in this work.Of 44, find that 13 (∼29%) show resolved morphology in optical filters.This is opposite to the general trend discovered above and implies that a fraction of these compact sources are likely experiencing a secondary burst, in a relatively small area, after the initial buildup of stellar structure.In particular, this is relevant to the stochastic nature of star formation in early galaxies, which may have a non-negligible effect on the observed UV magnitude measurements and consequently UV luminosity functions.

Summary
In this study, we have identified 341 galaxies at 5 < z < 14 in legacy fields of JWST and analyzed their rest-frame UV and optical sizes through JWST NIRCam images.The imaging data used here were collected from several public programs in Cycle 1, resulting in a combined effective area of 358 arcmin 2 .With a robust (8σ) selection of 341 galaxies, made possible by the unprecedented area coverage provided by JWST, we have conducted the first systematic exploration of the size-mass relation of galaxies in the first billion years.The key findings are as follows: 1.The slope of the size-mass relation was derived via linear regression analyses and found to be α ∼ 0.2, similar to those of star-forming galaxies at z < 3, but scaled down in size by 0.4 dex.The derived intercept was found to evolve with ∝(1 + z) −0.4 , a much slower evolution than those found in previous studies of rest-frame UV sizes of galaxies at lower redshifts.2. By using the results from our linear regression analysis, we identified 44 compact sources that are marginally resolved in NIRCam imaging.These compact sources account for ∼13% of the full sample presented here.3. We found that our sources overall have a high star formation surface density (1 M e yr −1 kpc −2 ), with the newly identified compact sources being as high as ∼300 M e yr −1 kpc −2 .We demonstrated that the absence of a clear declining trend indicates that the star formation efficiency may remain high even in the high-mass range; if the observed high efficiency remains similar in the following ∼0.5-1Gyr, some of our sources would evolve to ∼10 11 M e by z ∼ 5. 4. For 51 sources with available NIRSpec/MSA data, we investigated their ISM via the [O III]-to-Hβ line ratio.None of the compact sources are confidently classified as AGNs in the MEx diagram; however, the nature of the compact sources remains to be conclusively elucidated in a future study.A potential explanation for the observed high line ratios is high shock velocities, driven by intense star formation characterized by high Σ SFR . 5. We found that the sizes in the rest-frame and optical wavelengths are on average consistent.We attributed this to the short mass-doubling time (i.e., 1/sSFR < 100 Myr) of our sources, implying that they are actively building their structure coherently and are thus dominated by young stars.
With the unprecedented resolution and sensitivity provided by JWST, this work has demonstrated a comprehensive size analysis of galaxies in the first billion years.Of particular interest are the newly discovered compact populations.Specifically, the physical mechanisms that maintain the observed high star formation rate in such compact systems, which are likely under the strong influence of negative feedback, remain an open question.Recent JWST observations have identified a number of faint AGNs (Matthee et al. 2023;Onoue et al. 2023) and a complex of dusty AGN + young stellar populations (Akins et al. 2023;Furtak et al. 2023;Labbe et al. 2023), suggesting a greater prevalence of AGNs in high-z galaxies than previously thought.These emerging findings raise a caveat that the estimated star formation rate from our SED analysis may not accurately represent the intrinsic value, even though our spectroscopic analysis on the subsample did not find any immediate signatures of AGNs.Future spectroscopic follow-ups of the compact sources will provide further insights into their nature and their potential impacts in a broad cosmological context.NGDEEP has a comparable depth to JADESGDS, it has a relatively small area (a single NIRCam pointing).In addition, the NGDEEP field does not allow photometric selections at z < 9.7, the redshift range where the majority of the compact sources are identified.

Figure 1 .
Figure1.Radial flux profiles of webbpsf with various σ jitter values (dashed lines) generated for the JADESGDS F115W image, where the adopted profile (σ jitter = 0.022) is color-coded in cyan.The profile of an example bright star taken from the image is shown for comparison (red solid line).The inset shows the distribution of χ 2 /ν of the same star fitted by webbpsf for different σ jitter values.
is the best-fit total magnitude derived by galfit and m total is the total magnitude derived in Section 3.1, both measured in the rest-frame UV filter of interest for the target redshift range.Sources flagged in the galfit results are set C galfit = 1.The inferred physical properties are presented in Appendix B.

Figure 2 .
Figure 2. Sample distribution in the star formation rate-stellar mass plane (circles).Those with spectroscopic redshift (blue squares) and flagged as compact (magenta; Section 4.4) are marked accordingly.Those flagged in galfit results (Section 3.3) are shown by open symbols.The SFMS slope at z = 5 (black dashed lines; Speagle et al. 2014), the linear regression slope derived for our entire sample, ( ) * M M logSFR 0.81 log 10 0.31 8  = + (with the slope fixed to that of Speagle et al. 2014 at z = 5), and the derived range for the intrinsic scatter (0.37 dex; shaded regions) are shown.

Figure 3 .
Figure 3. Example of galfit fitting results in 2D images (left: original; middle: model; and right: residual).Top: PRIMERCOS-38203, one of the highest-redshift sources in our sample.Middle: JADESGDS-30934 (spectroscopically confirmed to z = 13.2;Curtis-Lake et al. 2023), an example of those classified as compact (Section 4.4).Bottom: F090W-dropout source that is flagged in the galfit fitting analysis.Flagged objects are not included in our statistical discussion.

Figure 4 .
Figure4.Distribution of our sample galaxies at 5 < z < 14 in the stellar mass-size plane (gray circles), in four redshift panels.The linear slope derived by regression analysis for the full sample (orange solid lines, with the shaded regions covering the 1.5R log e s range from the median slope; Table3) is shown.For comparison, the slope derived for late-type galaxies at z ∼ 2.7 (van der Wel et al. 2014; black dashed lines) is shown.Those classified as compact (Section 4.4) are shown by magenta symbols.Those with spectroscopic redshift are marked by open blue squares.The two horizontal hatched regions show the physical size of FWHM/2 for NIRCam (dark gray) and HST/WFC3-IR F160W (light gray) at the median redshift of the sources in each redshift window.

Figure 5 .
Figure 5.The same as in Figure 4, but for the redshift-corrected size ( ( ) R z 1 eff, maj

Figure 6 .
Figure6.Redshift evolution of rest-frame UV size of our galaxies (gray circles).The size of the individual galaxies shown here is corrected to the pivot mass (Section 4.2).The redshift evolution of the intercept, ( ) z 1 z µ + a , is shown (orange dashed line), along with those derived in two previous studies (the green line is forMosleh et al. 2012 and the black line is forvan der Wel et al. 2014).We note that the size trend of both previous studies is also scaled to the same pivot mass.Median sizes (large symbols) are derived in each redshift window.Those identified as compact (Section 4.4) are shown by magenta symbols.

Figure 7 .
Figure 7. Left: Star formation surface density (Σ SFR ) of our final sample (gray circles) as a function of redshift.The median values at each dropout redshift window are shown by larger symbols (squares).The upper limits for Σ SFR by spatial resolution limit, when assuming a star formation rate of 10 M e yr −1 , are shown for the corresponding JWST NIRCam filters (solid curved lines) and HST WFC3-IR F160W (black dotted line).The fit derived for LBGs in the HST data (Shibuya et al. 2015; cyan dashed line) is shown.Right: Σ SFR as a function of stellar mass.In the background, we show the density contour for the distribution of galaxies at 0.3 < z < 3.5, taken from the 3DHST catalog(Skelton et al. 2014;van der Wel et al. 2014).Those outside the lowest contour level are shown individually (crosses).Only ∼6% of the low-z sample is found at Σ SFR > 1 M e /yr/kpc 2 ; the majority of our sample is located above the same value.

Figure 8 .
Figure8.The same as Figure5, but the symbols are color-coded by the age ( t log ) derived in our SED analysis.The median age in each grid is shown by the larger symbols (squares), with the symbol size scaled by the number of galaxies in the grid.In the inset, the scatters in size measured at each mass bin are shown.The error of each scatter measurement is estimated by bootstrapping.Clearly, the radius and age of the stellar population increase, while the scatters in size decrease, with increasing stellar mass.

Figure 9 .
Figure 9. Left: distribution of 51 galaxies at 5 < z < 9.5 that have robust [O III] λ5007 /Hβ measurements by NIRSpec MSA, in the MEx diagram (circles).For those with Hβ detected at S/N <3, lower limits for the ratio are shown (triangles).Compact galaxies (as defined in Section 4.4; magenta) and two sources from Kocevski et al. (2023) and Larson et al. (2023; cyan stars) are marked separately.The solid lines are shown to indicate the regions dominated by star-forming galaxies (bottom left) and AGNs (top right) of z ∼ 2 (Juneau et al. 2014).Right: distribution of the normalized size ( R e R log e s D, an indicator that is used to define compactness in Section 4.4) as a function of the distance to the MEx border separation line.
Note.Limiting magnitudes measured in empty regions of the image with r = 0 16 apertures.a Mosaic images have been created using the first epoch data that are available as of 2023 June.
Note. "-d" represents dropout.The numbers of spectroscopically confirmed sources (see Section 3.2.3)are shown in brackets.

Table 3
Size-Mass Relations of Galaxies at 5 < z < 14: Linear Regression Best-fit Coefficients