#  Morphology of Galaxies in JWST Fields: Initial distribution and Evolution
of Galaxy Morphology

Jeong Hwan Lee Department of Astronomy and Atmospheric Sciences, Kyungpook
National University, Daegu 41566, Republic of Korea  The Center for High
Energy Physics, Kyungpook National University, Daegu 41566, Republic of Korea
Changbom Park Korea Institute for Advanced Study, 85 Hoegi-ro, Dongdaemun-gu,
Seoul 02455, Republic of Korea  Ho Seong Hwang Department of Physics and
Astronomy, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826,
Republic of Korea  SNU Astronomy Research Center, Seoul National University, 1
Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea  Minseong Kwon
Department of Physics and Astronomy, Seoul National University, 1 Gwanak-ro,
Gwanak-gu, Seoul 08826, Republic of Korea Ho Seong Hwang
<EMAIL_ADDRESS>
###### Abstract

A recent study from the Horizon Run (HR5) cosmological simulation has
predicted that galaxies with ${\rm log}~{}M_{\ast}/M_{\odot}\lesssim 10$ in
the cosmic morning ($10\gtrsim z\gtrsim 4$) dominantly have disk-like
morphology in the $\Lambda$CDM universe, which is driven by the tidal torque
in the initial matter fluctuations. For a direct comparison with observation
we identify a total of about $18,000$ James Webb Space Telescope (JWST)
galaxies with ${\rm log}~{}M_{\ast}/M_{\odot}>9$ at $z=0.6-8.0$ utilizing deep
JWST/NIRCam images of publicly released fields, including NEP-TDF, NGDEEP,
CEERS, COSMOS, UDS, and SMACS J0723$-$7327\. We estimate their stellar masses
and photometric redshifts with the dispersion of $\sigma_{\rm NMAD}=0.007$ and
outlier fraction of only about 5%. We classify galaxies into three
morphological types, ‘disks’, ‘spheroids’, and ‘irregulars’, applying the same
criteria used in the HR5 study. The morphological distribution of the JWST
galaxies shows that disk galaxies account for $60-70\%$ at all redshift
ranges. However, in the high-mass regime (${\rm
log}~{}M_{\ast}/M_{\odot}\gtrsim 11$), spheroidal morphology becomes the
dominant type. This implies that mass growth of galaxies is accompanied with
morphological transition from disks to spheroids. The fraction of irregulars
is about 20% or less at all mass and redshifts. All the trends in the
morphology distribution are consistently found in the six JWST fields. These
results are in close agreement with the results from the HR5 simulation,
particularly confirming the prevalence of disk galaxies at small masses in the
cosmic morning and noon.

Galaxy evolution (594) — Galaxy classification systems (582)

††software: Source Extractor (Bertin & Arnouts, 1996), GALFIT (Peng et al.,
2002, 2010), eazy-py (Brammer et al., 2008; Brammer, 2023b), WebbPSF (Perrin
et al., 2014), Numpy (Harris et al., 2020), Matplotlib (Hunter, 2007), Scipy
(Virtanen et al., 2020), Astropy (Astropy Collaboration et al., 2013, 2018,
2022)

##  1 Introduction

Since the Hubble’s tuning fork scheme for galaxy classification was introduced
(Hubble, 1926; Sandage, 1961), galaxy morphology has provided important
insights into the evolution of galaxies across cosmic time. Previous studies
have found that galaxy morphology is closely involved with the environment and
various intrinsic properties of galaxies such as luminosity, mass, and star
formation activity. For instance, the morphology-density and morphology-radius
relations (Dressler, 1980; Postman et al., 2005; Park et al., 2007; Park &
Hwang, 2009; Fasano et al., 2015) demonstrated that the proportion of early-
type galaxies (ellipticals and S0 galaxies) in galaxy clusters increases as
clustercentric distance decreases and galaxy number density increases. These
relationships are also dependent on the stellar mass of galaxies, showing
significantly high fractions of early-type galaxies ($\gtrsim 50\%$) for
galaxies more massive than ${\rm log}~{}M_{\ast}/M_{\odot}>11$ (morphology-
mass relation; Bamford et al., 2009; Vulcani et al., 2011; Calvi et al.,
2012). Furthermore, several studies using the Sloan Digital Sky Survey (SDSS)
data suggested that the morphology of neighboring galaxies also plays an
important role in the morphological distribution of galaxies in local universe
(Park et al., 2008; Park & Choi, 2009) and in the Great Observatories Origins
Deep Survey (GOODS) field at up to $z\sim 1$ (Hwang & Park, 2009). They showed
that the early-type probability of a galaxy increases when the galaxy is
closer to another early-type neighbor galaxy, whereas the probability
decreases when the neighbor is a late-type galaxy. In addition, Hwang & Park
(2009) showed that the morphology-density relation becomes much weaker at
$z\sim 1$ compared to the nearby universe, implying that the morphological
distribution of galaxies must depend on redshift.

Beyond $z>1$, high-resolution optical and near-infrared (NIR) images from the
Hubble Space Telescope (HST) have been used for investigating the morphology
of galaxies. Most studies utilizing HST data have consistently reported that
galaxies at $z\sim 2$ exhibit a high prevalence of disturbed or interacting
morphologies, which are categorized as peculiar types (Abraham & van den
Bergh, 2001; Conselice et al., 2005, 2008; Papovich et al., 2005; Huertas-
Company et al., 2009, 2016; Buitrago et al., 2013; Mortlock et al., 2013).

However, there have been two major problems in previous HST studies on high-
redshift galaxy morphology. First, there are notable discrepancies in the
specific trends of morphological distributions. Buitrago et al. (2013)
presented that the fraction of disk-like galaxies increases from approximately
$15\%$ at $z\sim 0$ to $\sim 80\%$ at $z>2$, whereas the fraction of spheroid-
like galaxies decreases from $\sim 85\%$ to $\sim 20\%$ over the same redshift
range, based on the data from the Palomar Observatory Wide-field
InfraRed/DEEP2 (POWIR/DEEP2) and GOODS NICMOS Survey (GNS) surveys. In
contrast, Mortlock et al. (2013) found a high fraction of spheroidal galaxies
($40\%$) and a negligible number of disk-like galaxies ($<10\%$) at $z>2$ in
the Ultra Deep Survey (UDS) region of the Cosmic Assembly Near-infrared Deep
Extragalactic Legacy Survey (CANDELS) field. These discrepancies could be
attributed to variations in the criteria for morphological classification and
sample selection. For instance, Buitrago et al. (2013) selected massive
galaxies (${\rm log}~{}M_{\ast}/M_{\odot}>11$) and classified their
morphologies mainly based on their Sérsic indices, while Mortlock et al.
(2013) only visually classified the morphology of galaxies with ${\rm
log}~{}M_{\ast}/M_{\odot}>10$. Second, HST optical and NIR data have
observational limitations for studying high-redshift galaxy morphology at
$z\gtrsim 3$. The longest wavelength coverage of the HST extends to $\sim
1.6{\rm\mu m}$ (WFC3-IR/F160W), restricting investigations of galaxy
morphology to $z\lesssim 2.5$ in the rest-frame optical wavelength range.
Furthermore, the spatial resolution of NIR (WFC3-IR) data is worse than that
of HST optical (ACS) data, showing a full width at half maximum (FWHM) twice
as broad ($\sim 0\farcs 2$) as that from the HST/ACS images. These limitations
pose challenges when studying the morphology of high-redshift galaxies.

Thanks to the unprecedented sensitivity and resolution of the James Webb Space
Telescope (JWST), there have been dramatic advances in observational studies
of high-redshift galaxy morphology. Recent JWST studies have utilized various
methods for classifying galaxy morphology, including visual inspections
(Ferreira et al., 2022, 2023; Jacobs et al., 2023; Kartaltepe et al., 2023),
supervised Convolutional Neural Networks (Huertas-Company et al., 2023), the
Morpheus deep-learning framework (Robertson et al., 2023), and unsupervised
machine learning (Tohill et al., 2023). Ferreira et al. (2022, 2023) conducted
visual classifications of galaxy morphology combined with nonparametric
measurements of concentration, asymmetry, and smoothness (CAS; Conselice,
2003), using JWST/NIRCam images of the SMACS 0723 and the Cosmic Evolution
Early Release Science Survey (CEERS) fields. They concluded that disk galaxies
are predominant at $z=3-6$, constituting $>40\%$ of galaxies with ${\rm
log}~{}M_{\ast}/M_{\odot}>9$. With a sample of 850 galaxies from the CEERS
field, Kartaltepe et al. (2023) also found a high fraction of disk galaxies,
comprising $60\%$ at $z\sim 3$ and $30\%$ at $z>6$, while the fractions of
spheroid and irregular galaxies are relatively constant at all redshifts.
Huertas-Company et al. (2023) and Tohill et al. (2023) investigated the
morphology of a relatively larger sample of galaxies at $z>2$ by employing
artificial intelligence techniques. They demonstrated that the fractions of
clumpy and irregular galaxies increase with increasing redshift, while the
fractions of spheroid-like galaxies decrease. Overall, recent JWST studies are
in agreement that the Hubble Sequence, encompassing diverse morphologies from
disk galaxies to spheroid galaxies, was already established in the early
universe at $z\sim 8-9$. This agreement results from the superior
observational performance of the JWST, which have improved the visibility of
regular morphological features such as disks, bulges, and spiral arms compared
to HST (Jacobs et al., 2023).

Figure 1:  Rest-frame wavelengths of JWST/NIRCam filters (F115W, F150W, F200W,
F277W, F356W, and F444W) as a function of redshift. The shaded gray region
represents the rest-frame optical wavelength range spanning $0.5-0.7~{}{\rm\mu
m}$.

In this work, we examine the rest-frame optical morphologies of a large sample
of $\sim 18,000$ galaxies at $z=0.6-8.0$ using data from six public fields of
JWST. We have improved the accuracy of photometric redshifts significantly,
which can potentially impact the study of redshift evolution of galaxy
morphology.  Figure 1 displays the rest-frame pivot wavelengths of JWST/NIRCam
filters (F115W, F150W, F200W, F277W, F356W, and F444W) as a function of
redshift. To probe the galaxy morphology, we choose the NIRCam filter to
observe the galaxies in rest-frame optical wavelengths ($0.5-0.7~{}{\rm\mu
m}$) at the corresponding redshifts. We then compare these galaxy morphologies
with the recent findings from Horizon Run 5 (HR5) cosmological simulation (Lee
et al., 2021). The HR5 results from Park et al. (2022) suggested that
approximately two-thirds of the galaxies with $M_{\ast}\leq
10^{10}~{}M_{\odot}$ in the cosmic morning ($10\gtrsim z\gtrsim 4$) exhibit
disk-like morphologies, implying the dominance of disks in the early universe.
The remaining galaxies are roughly equally divided into irregulars and
spheroids, with a slight prevalence of irregulars. Our objective is to verify
if these HR5 results are supported by observational evidence, using high-
redshift galaxies detected in JWST/NIRCam images when the morphology
classification criteria as close as those used in HR5 are applied. To achieve
this scientific goal, we selected galaxies with stellar mass of
$M_{\ast}>10^{9}~{}M_{\odot}$ across a wide range of redshifts ($z=0.6-8.0$).
In line with the HR5 study, we applied the same morphology classification
scheme as depicted in Figure 2 of Park et al. (2022).

Table 1: HST and JWST Images and Filters Used in This Study

Field | Telescope | Instrument | Filters | Coverage AreaaaAreas are measured based on the image observed by the instrument’s filter with the longest wavelength.
---|---|---|---|---
|  |  |  | (${\rm arcmin^{2}}$)
NEP-TDF | HST | ACS | F435W, F606W | 22.0
| JWST | NIRCam (SWC) | F090W, F115W, F150W, F200W | 10.5
|  | NIRCam (LWC) | F277W, F356W, F410M, F444W | 10.0
NGDEEP | HST | ACS | F435W, F606W, F775W, F814W, F850LP | 37.3
|  | WFC3-IR | F105W, F125W, F160W | 30.7
| JWST | NIRCam (SWC) | F115W, F150W, F200W | 11.2
|  | NIRCam (LWC) | F277W, F356W, F444W | 9.4
CEERS | HST | ACS | F435W, F606W, F814W | 201.3
|  | WFC3-IR | F125W, F140W, F160W | 167.4
| JWST | NIRCam (SWC) | F115W, F150W, F200W | 89.9
|  | NIRCam (LWC) | F277W, F356W, F410M, F444W | 91.7
COSMOS | HST | ACS | F435W, F475W, F606W, F814W | 245.6
|  | WFC3-IR | F125W, F140W, F160W | 241.3
| JWST | NIRCam (SWC) | F090W, F115W, F150W, F200W | 86.8
|  | NIRCam (LWC) | F277W, F356W, F410M, F444W | 88.9
UDS | HST | ACS | F435W, F606W, F814W | 202.8
|  | WFC3-IR | F125W, F140W, F160W | 186.6
| JWST | NIRCam (SWC) | F090W, F115W, F150W, F200W | 147.4
|  | NIRCam (LWC) | F277W, F356W, F410M, F444W | 150.8
SMACS0723 | HST | ACS | F606W, F814W | 12.3
|  | WFC3-IR | F105W, F125W, F140W, F160W | 5.5
| JWST | NIRCam (SWC) | F090W, F150W, F200W | 11.4
|  | NIRCam (LWC) | F277W, F356W, F444W | 11.0

††footnotetext: Note.

This paper is structured as follows. In Section 2.1, we explain the HST and
JWST images used in our analysis. Sections 2.2 and 2.3 describe how we
conducted multiwavelength photometry and estimated photometric redshifts and
stellar masses of galaxies through spectral energy distribution (SED) fitting.
Section 3 describes the application of the galaxy classification scheme, with
measurements of Sérsic indices and asymmetry. In Section 4, we illustrate our
major findings about the morphological fractions as a function of redshift and
stellar mass. In Section 5, we discuss the implications of our results,
including a comparison with the HR5 study. Finally, Section 6 summarizes our
key results. Throughout this paper, we adopt the cosmological parameters in
the Planck 2015 results (Planck Collaboration et al., 2016), which are also
employed in the HR5 simulation (Lee et al., 2021): $H_{0}=68.4~{}{\rm
km~{}s^{-1}~{}Mpc^{-1}}$, $\Omega_{M}=0.3$, $\Omega_{\Lambda}=0.7$,
$\Omega_{\rm b}=0.047$, and $\sigma_{8}=0.816$.

##  2 Data and Methods

###  2.1 HST and JWST/NIRCam Images

In this study, we made use of publicly available images of HST and JWST/NIRCam
from the Grizli Image Release
v6.0111https://s3.amazonaws.com/grizli-v2/JwstMosaics/v6/index.html
https://grizli.readthedocs.io/en/stable/grizli/image-release-v6.html, which
was reduced and processed by the Grizli pipeline (Brammer, 2023a). This
pipeline has been known to be effective to alleviate instrumental artifacts in
JWST countrate products (*_rate.fits), including vertical and horizontal
stripes from electronic read noise (commonly known as ‘$1/f$ noise’; Schlawin
et al., 2020), circular patterns by significant cosmic ray events
(‘snowballs’), and stray light features like ‘wisps’ and ‘claws’ (Rigby et
al., 2023). After correcting these artifacts, the pipeline aligned the HST and
JWST/NIRCam images using stars from GAIA DR3 catalog (Gaia Collaboration et
al., 2021), and combined those images. Detailed information about the image
reduction process was provided in previous JWST studies (Bezanson et al.,
2022; Bradley et al., 2023; Valentino et al., 2023). The final mosaic images
are sampled at pixel scales of $0\farcs 04~{}{\rm pixel^{-1}}$ for HST and
JWST/NIRCam long wavelength channels (F277W, F356W, and F444W) and $0\farcs
02~{}{\rm pixel^{-1}}$ for JWST/NIRCam short wavelength channels (F090W,
F115W, F150W, and F200W). For our analysis, we resampled the NIRCam short-
wavelength images to $0\farcs 04~{}{\rm pixel^{-1}}$ using the Python-based
tool astropy reproject222https://github.com/astropy/reproject.

We obtained the processed HST and JWST images of six public fields: the North
Ecliptic Pole Time-Domain Fields (NEP-TDF; Windhorst et al., 2022, 2023), the
Next Generation Deep Extragalactic Exploratory Public (NGDEEP) Survey
(Finkelstein et al., 2021; Bagley et al., 2023), the CEERS (Finkelstein et
al., 2017, 2023), the Cosmic Evolution Survey (COSMOS; Dunlop et al., 2021;
Kartaltepe et al., 2021; Casey et al., 2023), the UKIRT Infrared Deep Sky
Survey Ultra-deep Survey field (UDS; Dunlop et al., 2021), and the SMACS
J0723.3-7327 cluster (SMACS0723; Pontoppidan et al., 2022).  Table 1
summarizes the HST and JWST images used in our spectral energy distribution
(SED) analysis, excluding images from filters with either excessively narrow
coverage areas (e.g., HST F105W in most fields) or shallow depths (e.g., HST
F275W in CEERS and HST F435W in SMACS0723).

Figure 2 illustrates the surface brightness limit and coverage area of each
JWST field in various NIRCam filters. We measured the $1\sigma$ surface
brightness levels above the background within circular apertures with a radius
of $0\farcs 5$. For these measurements, we randomly selected 100 circular
apertures per field-of-view of $2\farcm 2\times 2\farcm 2$, corresponding to a
single JWST/NIRCam module. Overall, the JWST NIRCam images exhibit surface
brightness limits deeper than $\sim 25~{}{\rm mag~{}arcsec^{-2}}$ for the
short-wavelength filters and $\sim 26~{}{\rm mag~{}arcsec^{-2}}$ for the long-
wavelength filters, which even enables the detection of low surface brightness
galaxies (Ikeda et al., 2023). The NGDEEP ($\mu_{\rm lim}({\rm F444W})\sim
28.0~{}{\rm mag~{}arcsec^{-2}}$) and UDS (coverage of $\sim 150~{}{\rm
arcmin^{2}}$) are the deepest and widest fields in this study, respectively.

Figure 2:  Distributions of surface brightness limits and coverage areas for
JWST NIRCam filters (F115W, F150W, F200W, F277W, F356W, and F444W). Surface
brightness limits are measured within apertures with a radius of $0\farcs 5$.
Note that the SMACS0723 field has no JWST images from F115W.

###  2.2 Multiwavelength Photometry

For target selection and SED analysis for high-redshift galaxies, we performed
source detection and aperture photometry on HST and JWST images using Source
Extractor version 2.25.0333https://github.com/astromatic/sextractor
(SExtractor; Bertin & Arnouts, 1996). We employed the combined images of the
NIRCam long-wavelength filters (F277W+F356W+F444W) as the reference images for
performing forced photometry on both HST and JWST images. We utilized a two-
mode detection process with ‘cold’ and ‘hot’ modes as applied in Galametz et
al. (2013). The cold-mode configurations were adjusted to prevent
overdeblending of large and bright sources, whereas the hot-mode
configurations were set to detect faint objects close to bright sources.
Table 2 lists the input configuration parameters for the cold-mode and hot-
mode detections of SExtractor. In the case of SMACS0723 cluster, we adjusted
the configuration parameters slightly to obtain sharper background estimation
due to intense intracluster light. Subsequently, we combined the source
catalogs from the two-mode photometry, by retaining all sources from the cold
mode and supplementing them with hot-mode sources located outside the Kron
radii of cold-mode sources (Barden et al., 2012; Kartaltepe et al., 2023).
From the merged source catalogs, we initially selected galaxy candidates using
specific criteria based on colors and sizes as shown in  Figure 3. We applied
realistic color ranges ($-1.5<{\rm F200W-F277W}<1.5$ and $-1.5<{\rm
F277W-F356W}<1.5$) and $\texttt{FLAGS}\leq 4$ to remove artifacts. In
addition, we rejected point sources with compact half-light radii
($\texttt{FLUX\\_RADIUS}\leq 0\farcs 09$ or
$\texttt{FLUX\\_RADIUS}\leq-0.008\times(\texttt{MAG\\_AUTO}-25)+0\farcs 09$)
in the reference images. The numbers of all detected sources and initially
selected galaxy candidates are listed in the first and second columns of
Table 4.

Table 2: SExtractor Input Parameters Used for Multiwavelength Photometry

Parameter |  Cold-mode Configuration |  Hot-mode Configuration
---|---|---
DETECT_MINAREA |  20 |  20
DETECT_THRESH |  1.0 |  1.0
ANALYSIS_THRESH |  1.0 |  1.0
FILTER_NAME |  tophat_5.0_5x5.conv |  gauss_2.5_5x5.conv
DEBLEND_NTHRESH |  32 |  16
DEBLEND_MINCONT |  0.01 |  0.0001
SATUR_LEVEL |  100.0 |  100.0
BACK_SIZEaaParameter values within parentheses are applied for the SMACS0723 field. |  64 (32) |  16
BACK_FILTERSIZEaaParameter values within parentheses are applied for the SMACS0723 field. |  5 (3) |  3
BACK_PHOTOTHICKaaParameter values within parentheses are applied for the SMACS0723 field. |  48 (24) |  24 (12)

††footnotetext: Note.

Figure 3:  Left panel: Color-color diagram of the detected sources with
$\texttt{FLAGS}\leq 4$ across all JWST fields. Measurements of ${\rm
F200W}-{\rm F277W}$ and ${\rm F277W}-{\rm F356W}$ were conducted within the
Kron radii of the objects. The blue dashed box denotes the selection criteria
of $-1.5<{\rm F200W}-{\rm F277W}<1.5$ and $-1.5<{\rm F277W}-{\rm F356W}<1.5$
to discard artifacts. Right panel: Size-magnitude diagram of the sources, with
half-light radii (FLUX_RADIUS) and magnitudes measured in the reference image
(${\rm F277W}+{\rm F356W}+{\rm F444W}$). The blue dashed line represents the
boundary between extended sources and point sources. Point source candidates
located below this boundary were excluded from our sample for the analysis.

###  2.3 Photometric Redshifts and Stellar Masses

####  2.3.1 SED Fitting Based on eazy-py

Prior to conducting SED fitting for galaxy candidates, we carried out an
aperture correction procedure to determine the total flux of each source,
following the method described in Valentino et al. (2023). We measured the
enclosed flux within a circular aperture with a radius of $0\farcs 5$ and
accounted for additional flux beyond the $0\farcs 5$ aperture within the Kron
radius (Kron, 1980). These aperture-corrected flux values were computed from
the reference images and were subsequently applied to all HST and JWST
filters. We employed these aperture-corrected total fluxes for the SED
analysis described in Section 2.3. To correct galactic extinctions in all HST
and JWST bands, we made use of dustmaps (Green, 2018), with adopting a
standard $R_{V}=3.1$ and dust reddening maps provided in Schlegel et al.
(1998).

Photometric redshifts and stellar masses of the galaxy candidates are
estimated through the SED fitting process using the eazy-
py444https://github.com/gbrammer/eazy-py package (Brammer et al., 2008;
Brammer, 2023b). Multiwavelength fluxes measured in HST and JWST bands (see
Table 1) were used for fitting the eazy-py model. We set the flux
uncertainties to a minimum error floor of $10\%$ to allow some flexibility in
the SED fitting. We employed 14 spectral templates from the sfhz
model555https://github.com/gbrammer/eazy-photoz/tree/master/templates/sfhz,
including 13 templates in the corr_sfhz_13 subsets and an additional best-fit
template from a strong emission-line galaxy observed in the SMACS0723 field
(Carnall et al., 2023). Photometric redshifts were determined by identifying
the redshift value corresponding to the minimum chi-sqaure value within the
range from $z=0.05$ to $z=12.0$, with an interval of ${\rm
log}~{}(1+z)=0.005$. An iterative correction process for flux zeropoints was
performed over 10 steps to optimize redshift values. No priors were applied on
apparent magnitude and ultraviolet (UV) slope during this procedure.

Figure 4:  Examples of the SED fitting results for six galaxies in the CEERS
field, ordered by their chi-square values of the best-fit eazy-py models. In
this figure, there are three panels for each galaxy. The upper panel in each
set displays the SED of the best-fit model (gray solid line) and photometric
data from HST (orange triangles) and JWST (magenta circles) observations. The
lower panel shows the chi-square values as a function of redshift. The
spectroscopic redshift and the photometric redshift are depicted as red and
blue dashed lines, respectively. A small panel on the right side displays the
observed galaxy in the reference image (F277W+F356W+F444W). Based on a
criterion of $\chi^{2}=20$, four galaxies in blue boxes are considered to show
secure SED fitting results, while two galaxies in red boxes are thought to
show unreliable results.

After completing the SED fitting process, we evaluated the reliability of the
best-fit SED models by utilizing the chi-square values ($\chi^{2}_{\rm SED}$)
from eazy-py and the number of filters ($N_{\rm filt}$) used for the SED
analysis. Through visual inspection of the SED fitting results, we employed a
criterion of $\chi^{2}_{\rm SED}<20$ and $N_{\rm filt}\geq 6$ to select well-
fitted galaxy candidates. It is worth noting that our chi-square criterion is
more rigorous compared to that applied in Valentino et al. (2023), which used
$\chi^{2}_{\rm SED}/N_{\rm filt}\leq 8$ and $N_{\rm filt}\geq 6$.  Figure 4
illustrates examples of SED fitting results for six galaxies located in the
CEERS field. Among these, four galaxies with $\chi^{2}_{\rm SED}<20$ exhibit
well-matched SEDs with their photometric data, except for the low-signal F606W
data for the galaxy of ID 6910. In contrast, the best-fit SEDs for two
galaxies with $\chi^{2}_{\rm SED}\geq 20$ do not align with the photometric
data. These poor results might be associated with contamination from a
neighboring bright galaxy (ID 80324) and the faintness of the galaxy itself
(ID 41179), as seen in visual representations of the galaxies in the figure.
For further analysis, we decided to exclude samples with $\chi^{2}_{\rm
SED}\geq 20$ or $N_{\rm filt}<6$ from the eazy-py SED fitting.

Figure 5:  Comparison of photometric redshifts ($z_{\rm phot}$) derived from
our eazy-py model in this study with spectroscopic redshifts ($z_{\rm spec}$)
obtained from the NED and the literature. The left panel displays the
comparison results for all sources, and the right panel displays the results
for sources selected with ${\rm log}~{}M_{\ast}/M_{\odot}>9$, $\chi^{2}_{\rm
SED}<20$, and $N_{\rm filt}\geq 6$. Gray dashed lines indicate the criterion
of $|z_{\rm spec}-z_{\rm phot}|/(1+z_{\rm spec})=0.15$, which is used to
identify outliers. This figure provides the total numbers of sources, the
numbers of outliers, the outlier fractions, and the $\sigma_{\rm NMAD}$ values
as text labels.

####  2.3.2 Comparisons of Photometric and Spectroscopic Redshifts

In this section, we compare photometric redshifts with spectroscopic redshifts
to evaluate the reliability of our photometric redshift estimations. The
spectroscopic redshift data for galaxies within JWST fields were obtained from
the NASA/IPAC Extragalactic Database (NED) and the available literature. In
the SMACS0723 field, the spectroscopic redshifts were sourced from Noirot et
al. (2023), which provided the first redshift catalog for this field utilizing
JWST/NIRISS grism, JWST/NIRSpec, and VLT/MUSE data. For the other JWST fields,
we obtained spectroscopic redshifts through a cone search in
NED666https://ned.ipac.caltech.edu/conesearch. These spectroscopic sources
were spatially matched with our photometric sources with ${\rm
F200W}<27~{}{\rm mag}$, allowing for a tolerance of $1\farcs 5$. This
tolerance was determined to consider potential offsets in the world coordinate
system (WCS) between JWST images and previous data. In the NEP-TDF field,
there was only one object with measured spectroscopic redshift (WISEA
J172302.06+654802.9; $z_{\rm spec}=0.179$), which is statistically
insufficient for redshift comparisons.

Figure 5 displays the comparisons of photometric and spectroscopic redshifts
from five JWST fields except for NEP-TDF. We plot the results from all matched
sources in the left panel and those from the sources with stellar mass higher
than ${\rm log}~{}M_{\ast}/M_{\odot}>9$, low chi-square values ($\chi^{2}_{\rm
SED}<20$), and a minimum six available filters ($N_{\rm filt}\geq 6$) in the
right panel. Outliers in this comparison were identified using the criterion
of $|z_{\rm spec}-z_{\rm phot}|/(1+z_{\rm spec})>0.15$ (Weaver et al., 2023),
which are located outside the gray dashed lines in the figure. To assess the
accuracy of photometric redshifts, we also computed the normalized median
absolute deviation ($\sigma_{\rm NMAD}$) of the redshift differences ($\Delta
z=z_{\rm phot}-z_{\rm spec}$) as defined in Equation 7 of Brammer et al.
(2008). For all matched sources, the photometric redshifts of galaxy
candidates show a good agreement with their spectroscopic redshifts in the
redshift range from $z\sim 0.1$ to $z\sim 10$, with an outlier fraction below
$8\%$ and $\sigma_{\rm NMAD}=0.008$. When the stricter criterion of ${\rm
log}~{}M_{\ast}/M_{\odot}>9$, $\chi^{2}_{\rm SED}<20$, and $N_{\rm filt}\geq
6$ is applied, the results become more reliable, yielding a lower outlier
fraction of $\sim 5\%$. The performance of photometric redshift estimation
slightly varies with fields, as detailed in  Table 3. The outlier fractions
and $\sigma_{\rm NMAD}$ values from our analysis are lower than those in
previous JWST studies, which reported an outlier fraction of $10.8\%$ and
$\sigma_{\rm NMAD}=0.03$ in the UNCOVER survey field (Weaver et al., 2023),
and $\sigma_{\rm NMAD}=0.018$ in the CEERS field (Valentino et al., 2023).
This implies that our sample is more rigorously selected based on the SED
fitting results than in previous studies.

Table 3: Comparison Results of $z_{\rm phot}$ versus $z_{\rm spec}$ in Each
JWST Fielda, ba, bfootnotemark:

Field | Number of Matched Sources | Number of Outliers | Outlier Fraction (%) | $\sigma_{\rm NMAD}$
---|---|---|---|---
| All $z$ | Selected $z$ | All $z$ | Selected $z$ | All $z$ | Selected $z$ |
NGDEEP | 51 (37) | 31 (26) | 6 (3) | 3 (3) | 11.8 (8.1) | 9.7 (11.5) | 0.010
CEERS | 375 (296) | 224 (204) | 28 (14) | 20 (12) | 7.5 (4.7) | 12.5 (6.9) | 0.007
COSMOS | 638 (541) | 472 (420) | 49 (24) | 44 (23) | 7.7 (4.4) | 9.3 (5.4) | 0.009
UDS | 192 (148) | 160 (129) | 14 (5) | 9 (5) | 7.3 (3.4) | 5.6 (3.9) | 0.006
SMACS0723 | 185 (159) | 72 (63) | 16 (8) | 7 (3) | 8.6 (5.0) | 9.7 (4.8) | 0.011
Total | 1,441 (1,181) | 998 (868) | 113 (54) | 86 (46) | 7.8 (4.6) | 8.6 (5.3) | 0.008

††footnotetext: Notes.

aafootnotetext: Values in parentheses correspond to results with ${\rm
log}~{}M_{\ast}/M_{\odot}>9$, $\chi^{2}<20$, and $N_{\rm filt}\geq 6$.

bbfootnotetext: The selected range of redshift is $z=1.2-8.0$ for SMACS0723
and $z=0.6-8.0$ for the other fields.

However, it is notable that there are several sources at $z_{\rm spec}\lesssim
1$ showing overestimated photometric redshifts with $z_{\rm phot}\gtrsim 1$.
This discrepancy can be attributed to contamination from adjacent light
sources, such as blending with other galaxy light, bright background features,
or bleeding light from saturated stars. These overestimated cases could be
low-redshift contaminants in our analysis. To minimize these contaminants, we
manually excluded a few sources close to saturated stars and bright galaxies
for all JWST fields through visual inspections.

Figure 6:  Comparisons of $z_{\rm phot}$ derived from the EAZY-PY model (left
panel; this study) and the RELICS catalog (right panel; Coe et al., 2019)
adopted by Ferreira et al. (2022) with $z_{\rm spec}$ (Noirot et al., 2023),
for the sources detected in the SMACS0723 field. Gray and purple lines and
texts are the same as  Figure 5.

In addition, we also compared the accuracy of photometric redshifts of
galaxies in the SMACS0723 field with those listed in the RELICS catalog (Coe
et al., 2019). The RELICS study derived photometric redshifts using the
Bayesian photoz code (BPz; Benítez, 2000), with HST photometric data from
optical to NIR wavelengths. These photometric redshifts were adopted by
Ferreira et al. (2022) for investigating the morphologies of high-redshift
galaxies.  Figure 6 represents the comparisons of $z_{\rm phot}$ vs. $z_{\rm
spec}$ both for this study (left panel) and the RELICS catalog (right panel).
The RELICS sources have 180 matched objects with the spectroscopic sample from
Noirot et al. (2023), which is similar to this study. In comparison to the
RELICS catalog, our photometric redshift measurements demonstrate a superior
match with spectroscopic redshifts, showing a lower outlier fraction and
$\sigma_{\rm NMAD}$. On the other hand, the RELICS catalog tends to
underestimate redshifts of some galaxies at $z_{\rm spec}>1$, which might be
due to observational limitations in HST NIR data capturing the redshifted
Balmer breaks of these galaxies. These suggest that some high-redshift galaxy
candidates in the SMACS0723 field might have been missed in Ferreira et al.
(2022) because of underestimation of photometric redshifts in the RELICS data.
Furthermore, our analysis demonstrates that the inclusion of JWST NIRCam
photometric data improves the accuracy of photometric redshifts when compared
to the results obtained using HST-only data.

Figure 7:  Distribution of stellar masses of extended sources in all JWST
field as a function of redshift. Stellar masses and redshifts are derived
using the EAZY-PY model in this study. The blue dashed lines represent the
selection criteria of ${\rm log}~{}M_{\ast}/M_{\odot}>9$ and $z=0.6-8.0$,
which are used to select the galaxy candidate for morphology classification.

####  2.3.3 Sample Selection Based on Redshifts and Stellar Masses

This section describes our sample selection from the results of SED fitting.
Initially, we established a robust sample by excluding the sources with high
chi-square values ($\chi^{2}_{\rm SED}\geq 20$), smaller number of filters
($N_{\rm filt}<6$), and implausible output parameters such as negative values
for photometric redshift, chi-square, stellar mass, and rest-frame fluxes. In
addition, we rejected faint objects with ${\rm F200W}>28~{}{\rm mag}$ from our
sample due to their low signal-to-noise ratios, which could be inappropriate
for the investigation of morphologies (Huertas-Company et al., 2023;
Kartaltepe et al., 2023). The numbers of these robust samples from the SED
fitting are listed in the third column in  Table 4. We then selected the high-
redshift galaxy samples for the morphological analysis based on redshifts and
stellar masses.  Figure 7 illustrates the stellar mass distribution for the
selected objects across all JWST fields, as a function of redshift. We
utilized spectroscopic redshifts when available, otherwise we employed
photometric redshifts obtained through SED fitting. Our final sample for
studying morphologies was chosen with a redshift range from $z=0.6$ to $z=8.0$
and a stellar mass range of ${\rm log}~{}M_{\ast}/M_{\odot}>9$, as represented
in the figure. For the galaxies in the SMACS0723 field, we carefully selected
a sample at the redshift range of $z=1.2-8.0$ to minimize any potential
contamination from the foreground SMACS0723 cluster. These selection criteria
ensure that our sample is not significantly affected by completeness issues of
observational data in the selected redshift range. The stellar mass range
chosen in this study was commonly employed in most JWST studies focusing on
galaxy morphology (Ferreira et al., 2023; Huertas-Company et al., 2023;
Kartaltepe et al., 2023). The specific number of objects selected based on
redshifts and stellar masses for each field is detailed in the fourth column
in  Table 4.

Table 4: Numbers of Objects in Each Field for All Selection Steps

Field |  Number of objects for each step
---|---
|  (1) |  (2) |  (3) |  (4) |  (5)
NEP-TDF |  12,222 |  4,929 |  2,314 |  799 |  755
NGDEEP |  17,184 |  6,276 |  1,971 |  551 |  501
CEERS |  86,484 |  44,956 |  13,578 |  5,159 |  4,991
COSMOS |  65,430 |  35,620 |  9,818 |  4,118 |  4,031
UDS |  89,178 |  54,101 |  21,174 |  7,607 |  7,408
SMACS0723 |  11,930 |  6,532 |  2,485 |  540 |  364
Total |  282,428 |  152,414 |  51,340 |  18,774 |  18,050

††footnotetext: Notes. Columns represent the numbers of the following source
categories: (1) The total number of sources detected by SExtractor; (2)
Extended source candidates chosen for the eazy-py SED fitting, meeting the
criteria based on colors and sizes as described in Section 2.2; (3) Robust
samples obtained from the eazy-py SED fitting, based on criteria such as
$0<\chi^{2}_{\rm SED}<20$, $N_{\rm filt}\geq 6$, ${\rm F200W}<28~{}{\rm mag}$,
and positive values for photometric redshifts, stellar masses, and rest-frame
fluxes. These sources are depicted in  Figure 7; (4) Galaxy candidates with
${\rm log}~{}M_{\ast}/M_{\odot}>9$ at $z=0.6-8.0$; (5) The final sample
selected through visual inspections and excluding non-converged sources with
$n<0.2$, $n>10$, $n\leq\delta n$, $|x_{0}^{\rm GALFIT}-x_{0}^{\rm
SExtractor}|>10~{}{\rm pixels}$, or $|y_{0}^{\rm GALFIT}-y_{0}^{\rm
SExtractor}|>10~{}{\rm pixels}$. We examine the morphologies of this final
sample using GALFIT, as described in Section 4.

##  3 Morphology Classification

###  3.1 Morphology Classification Scheme in the HR5 Study

For fair comparisons of the results on galaxy morphology among different
studies, the classification scheme should be as close as possible. To ensure
consistency with the HR5 study, we adopted the same morphology classification
scheme described in Park et al. (2022). In their work, approximately $33,000$
simulated galaxies with $M_{\ast}>2\times 10^{9}~{}M_{\odot}$ were categorized
into three distinct types: ‘disks’, ‘spheroids’, and ‘irregulars’. As depicted
in Figure 2 of Park et al. (2022), their classification scheme relies on two
key parameters: the Sérsic index ($n$) and the asymmetry factor ($A$). These
parameters were derived from the spatial distribution of stellar particles
belonging to each galaxy.

Initially, symmetric galaxies were categorized with $A<0.4$, while asymmetric
galaxies were characterized by $A>0.4$. This criterion is similar to
thresholds used in previous observational studies (Conselice, 2003; Ferreira
et al., 2023; Kartaltepe et al., 2023), which differentiated undisturbed
galaxies with $A<0.35$ and merging galaxies with $A>0.35$. However, the
asymmetry measurement in Park et al. (2022) was based on the three-dimensional
distribution of stellar particles, which could lead to a systematic difference
from the asymmetry calculated using two-dimensional observational data. In
this study, we adopted a criterion of $A=0.32$ derived from the direct
comparison between three- and two-dimensional asymmetry parameters (refer to
Appendix B for details).

In the HR5 study, the radial stellar mass distribution within galaxies was fit
by the Sérsic model. Symmetric galaxies were further divided into disks disks
($A<0.4$ and $n<1.5$) and spheroids ($A<0.4$ and $n>1.5$) based on the results
of Sérsic profile fitting. The choice of the Sérsic index threshold of $n=1.5$
was based on the fact that it effectively distinguishes between late-type and
early-type galaxies, as elaborated in the Appendix of Park et al. (2022). It
is very important to note that the criterion $n=1.5$ results in the most
accurate classification when the central part of galaxies is excluded in the
Sérsic profile fitting as this can remove the galaxy bulge component and any
observational artifact/AGN near the center. When the central part is included
in the fitting, disk galaxies with a bright central bulge or core star-burst
will become difficult to be distinguished from spheroidal ones.

Irregular galaxies were categorized to include asymmetric galaxies ($A>0.4$)
or galaxies where the Sérsic fitting failed. In the subsequent sections, we
explain how we applied this classification scheme to categorize galaxy
morphologies using JWST observational data.

###  3.2 Measuring the Asymmetry

The rotational asymmetry index has commonly used for quantifying the degree of
asymmetry in galaxy light, as introduced by Abraham et al. (1996). This index
is determined by calculating the sum of differences between the original image
and its corresponding $180^{\circ}$-rotated image, relative to the sum of the
original image. In Park et al. (2022), the asymmetry of simulated galaxies was
measured using the three-dimensional stellar mass density within pixels
containing at least 10 stellar particles, excluding the central region within
a radius of $0.8~{}{\rm kpc}$ in proper distance.

In this work, we followed a similar methodology to Park et al. (2022) with the
additional correction for observational background noise, as recommended in
Conselice et al. (2000) and Lotz et al. (2004). The asymmetry index ($A$) was
computed using the following formula,

$A=\frac{\Sigma_{i,j\in S}~{}|I_{0}(i,j)-I_{180}(i,j)|}{2\times\Sigma_{i,j\in
S}~{}I_{0}(i,j)}-A_{\rm bkg},$ (1)

where $I_{0}$ and $I_{180}$ are the pixel values in the original and
$180^{\circ}$-rotated images, and $A_{\rm bkg}$ is the asymmetry index for the
background region. The indices of $i$ and $j$ denote the pixels within the
selected region ($S$) used to measure asymmetry.

To apply the formula, we followed these steps. Initially, we created cutout
images for all selected galaxies with a size of $10^{\prime\prime}\times
10^{\prime\prime}$, which is sufficient to cover the total fluxes of the
galaxies. The filter selection for these cutout images was tailored to the
redshift of each galaxy, as illustrated in  Figure 1. Then, we ran SExtractor
twice on the cutout images, with a reduced background mesh size of 32, and
keeping other input configurations consistent with those listed in  Table 2.
The segmentation map from the initial SExtractor run was used to mask any
neighboring light sources adjacent to the target galaxy. This masking process
was implemented by replacing the pixel values of the detected contaminants
with the background value of the image, creating a masked cutout image. In the
second SExtractor run with the masked cutout image, we were able to more
precisely measure the Petrosian radius ($r_{p};$ Petrosian, 1976) of each
galaxy. We also generated a local background image around the galaxy within
the masked image, by utilizing a two-dimensional linear function to fit any
possible background gradients. Although the obtained JWST images were already
background-subtracted, these cutout images might still contain subtle
background gradients. This background gradient model was used in the
subsequent GALFIT analysis described in Section 3.3.

To calculate the asymmetry parameter, we selected a specific region for the
detected source ($S$) within an elliptical aperture with a semi-major axis of
$1.5\times r_{p}$ from the galaxy center. The radius of $1.5\times r_{p}$ is
known to encompass nearly the total flux of a galaxy (Conselice et al., 2000;
Lotz et al., 2004). The inner circular region with a radius of $0.8~{}{\rm
kpc}$ was excluded for this calculation, as done in Park et al. (2022). The
pixel values within the region $S$ were summed without any smoothing
procedures. To correct for the background noise, we designated the background
region as the area outside a radius of $1.5\times r_{p}$ from the center of
the galaxy. For determining $A_{\rm bkg}$, we adopted the minimum asymmetry
value from the background regions in the $90^{\circ}$-, $180^{\circ}$-, and
$270^{\circ}$-rotated images. This minimum asymmetry value was then normalized
by the area of the object region ($S$) when applying Equation 1.

###  3.3 Surface Photometry with GALFIT

We utilized GALFIT (Peng et al., 2002, 2010) to analyze the morphology of
galaxies based on their light distribution. We used the masked cutout image of
each galaxy, as mentioned in Section 3.2, as the input image for running
GALFIT. We applied a single Sérsic profile with free parameters for the
central coordinates, total magnitude, effective radius, and Sérsic index. The
axis ratio and position angle were fixed to the output values derived from the
second SExtractor run described in Section 3.2, in order to enhance the
stability of GALFIT solutions and speed up the fitting procedure. In addition,
we applied a two-dimensional sky component with a free sky background value
and fixed sky gradients along the x and y axes.

It is important to note certain limitations in our choice of a single Sérsic
profile for GALFIT in this study. Galaxies are known to exhibit multiple
morphological components such as disks, bars, spiral arms, bulges, and
irregular features resulting from merger or tidal processes. In particular,
around $10\%$ of galaxies have both bulge components ($n\sim 4$) at their
centers and exponential disk components ($n\sim 1$) in their outskirts
simultaneously (Huertas-Company et al., 2023; Kartaltepe et al., 2023). This
inherent complexity in galaxy morphology has led previous studies on high-
redshift galaxy morphology to rely on visual classification and use parametric
measurements from GALFIT as a supplement. Ideally, the inclusion of multiple
Sérsic profiles in the GALFIT configurations would be more beneficial for
tracing the detailed morphological structures of real galaxies. However,
multiple Sérsic fitting does not work effectively for most high-redshift
galaxies with small apparent size ($R_{\rm eff}\lesssim 0\farcs 2$), and it is
very challenging to apply to a large sample of galaxies. Therefore, this work
opted to use a single Sérsic profile for the consistency in methodology with
the HR5 study and the convenience of analyzing a large sample of high-redshift
galaxies.

In our GALFIT configurations, we employed point spread function (PSF) images
generated from WebbPSF (Perrin et al., 2014) to account for the PSF effect. We
resampled the PSF images of NIRCam bands to the pixel scale of $0\farcs
04~{}{\rm pixel^{-1}}$. The size of the PSF convolution box was set to be
larger than the total area of the input images. To minimize the cases of
obtaining unreasonable results, we imposed constraints on the Sérsic index,
ranging from 0.2 to 10, and the coordinate offsets to $0\farcs 1$ (2.5 pixels)
from the galaxy center.

After executing GALFIT, we systematically rejected sources with unreliable
solutions, such as large deviations in central coordinates ($|x_{0}^{\rm
GALFIT}-x_{0}^{\rm SExtractor}|>10~{}{\rm pixels}$, $|y_{0}^{\rm
GALFIT}-y_{0}^{\rm SExtractor}|>10~{}{\rm pixels}$), or the uncertainty of the
Sérsic index ($\delta n$) exceeding the value of the Sérsic index itself
($n\leq\delta n$). The final sample size with this rejection is detailed in
the fifth column of  Table 4. We additionally evaluated the quality of the
GALFIT fitting by computing the residual flux fraction ($RFF$; Hoyos et al.,
2011) for each galaxy. The $RFF$ measures the signal in a residual image
relative to the sum of the original image. As in Hoyos et al. (2012), the
$RFF$ is defined as,

$RFF=\frac{\Sigma_{i,j\in S}~{}|I(i,j)-I(i,j)^{\rm
model}|-0.8\times\Sigma_{i,j\in S}~{}\sigma_{\rm B}(i,j)}{\Sigma_{i,j\in
S}~{}I(i,j)}$ (2)

where $|I(i,j)-I(i,j)^{\rm model}|$ is the residual signal between the input
image and the GALFIT model image. $\sigma_{\rm B}$ denotes the background
fluctuation, and $S$ is the region selected for this calculation. In this
study, we selected $S$ as an elliptical region within a semi-major axis of
$1.5\times r_{p}$, excluding the central region within a radius of 0.8 kpc. In
terms of the background fluctuation, we adopted the background sigma for the
total area of the masked cutout image, as assumed in Margalef-Bentabol et al.
(2016). The equation used for calculating the background fluctuation is as
follows,

$\Sigma_{i,j\in S}~{}\sigma_{B}(i,j)=N_{S}\times\langle\sigma_{\rm B}\rangle,$
(3)

where $N_{S}$ is the number of pixels belonging to the area used for $RFF$
calculation, and $\langle\sigma_{B}\rangle$ is the mean background sigma in
the image. We determined $\langle\sigma_{B}\rangle$ to be the clipped sigma
value from the region outside the Petrosian aperture ($S$).

We utilized the $RFF$ as a proxy to assess the success of a single Sérsic
fitting. The specific $RFF$ criterion has not been agreed well because it can
vary depending on the quality of observational data and the scientific
purposes of using the $RFF$. In a recent JWST study by Ormerod et al. (2023)
and Ward et al. (2023), the $RFF$ was employed as a criterion for excluding
poorly-fit objects with $RFF>0.5$. This threshold was found to be more
generous than previous HST studies, which suggested that objects with
$RFF>0.11$ are needed to fit with additional Sérsic components (Hoyos et al.,
2011), or merger candidates could be effectively found with $RFF>0.2$
(symmetric cases) or $RFF>0.1$ (asymmetric cases). In this work, we adopted
the criterion of $RFF>0.5$ to identify obviously poorly-fit objects from the
GALFIT analysis along with avoiding any misclassification of multi-Sérsic
objects as poorly-fit ones.

###  3.4 Parametric Morphological Classification in This Study

Figure 8:  Thumbnail images ($6^{\prime\prime}\times 6^{\prime\prime}$) of
high-redshift galaxies selected in this study. The examples include six
galaxies, with two galaxies in each morphological category: disks (upper two
objects), spheroids (middle two objects), and irregulars (lower two objects).
Each galaxy example is represented by five panels. The leftmost panels display
the masked cutout images used as input for GALFIT. The subsequent two panels
display the model images and residual images obtained from GALFIT. In the
fourth panels, the mask images show the distribution of masked pixels in the
gray region, with the black region representing the $1\sigma$ boundary of
light from objects recorded in the segmentation map. Purple ellipses represent
apertures with a semi-major axis of $1.5\times r_{p}$ of each galaxy. The
fifth panels show the radial light profiles measured from the galaxy center to
the image edge. Original light profiles are depicted by blue curves, while the
model light profiles from GALFIT are depicted by red curves. Key parameters
for morphology classification are denoted by blue texts in this figure.

Utilizing the results obtained from asymmetry measurements and GALFIT fitting,
we applied the following parametric criteria for morphology classification:

(1) Disk-type galaxies: $n<1.5$, $A<0.32$, and $RFF<0.5$,

(2) Spheroid-type galaxies: $n>1.5$, $A<0.32$, and $RFF<0.5$,

(3) Irregular-type galaxies: $A>0.32$ or $RFF>0.5$.
Symmetric and well-fitted galaxies with $A<0.32$ and $RFF<0.5$ were further
categorized into disk-type ($n<1.5$) and spheroid-type galaxies ($n>1.5$)
based on their Sérsic indices. Irregular-type galaxies in our analysis
encompass both asymmetric galaxies with $A>0.32$ and poorly-fit galaxies with
$RFF>0.5$. This parametric classification scheme is consistent with the
criteria applied in Park et al. (2022).

Figure 8 illustrates the examples of our morphological classification. Disk
galaxies in the top two panels show nearly exponential light profiles ($n\sim
1$), occasionally accompanied by spiral arm structures (ID 22014). Spheroid
galaxies in the middle two panels are characterized by bulge-like structures
with higher Sérsic indices ($n\gtrsim 2$) compared to disk galaxies. In the
bottom two panels, irregular galaxies exhibit more remarkable substructures in
their residual images compared to disks and spheroids. Although there could be
a potential problem of misjudgement near the sharp classifying boundary of
parameters, our parametric classification scheme generally aligns with the
common visual classification of galaxy morphology for the majority of galaxies
in our extensive sample.

##  4 Results

###  4.1 Distributions of Asymmetry, $RFF$, and Sérsic indices

Figure 9:  Histograms depicting the asymmetry of galaxies in all JWST fields.
The panels are organized based on the chosen JWST/NIRCam filters corresponding
to the respective redshift ranges. The gray dashed lines represent our
criterion to distinguish symmetric galaxies ($A<0.32$) and asymmetric
($A>0.32$) galaxies. The black lines mark the median values of asymmetry in
each redshift bin.

In this section, we present the distributions of key parameters for morphology
classification: asymmetry, $RFF$, and Sérsic indices.  Figure 9 shows the
distributions of asymmetry values of our sample galaxies. The asymmetry
distributions exhibit a strong concentration around $A\sim 0.18$ across all
redshift ranges. As also shown in  Table 5, the median asymmetry values of
galaxies in each JWST field are also around $A\sim 0.18$, indicating no
systematic biases in asymmetry measurements depending on filters and fields.
Figure 9 also shows elongated tails towards higher asymmetry values in all
redshift ranges. These tails are indicative of irregular galaxies with mergers
and remarkable substructures. As mentioned in Section 3.1, we classified
galaxies as asymmetric when $A>0.32$, which is consistent with the criterion
of three-dimensional asymmetry with $A>0.4$ in Park et al. (2022).

Figure 10:  Distributions of $RFF$ values as a function of asymmetry. The panels are ordered as in  Figure 9. The gray lines denote the selection criteria of $RFF>0.5$ and $A>0.32$ to identify irregular galaxies in this study. The gray cross symbols mark the median values of $RFF$ and asymmetry in each redshift bin.  Table 5: Statistics of asymmetry, $RFF$, and Sérsic index of Galaxies in Each JWST Field Field | Asymmetry |  $RFF$ | Sérsic index
---|---|---|---
| Median |  $\sigma_{\rm MAD}$ | Median |  $\sigma_{\rm MAD}$ | Median | $\sigma_{\rm MAD}$
NEP-TDF | 0.18 | 0.11 | 0.14 | 0.11 | 1.06 | 0.39
NGDEEP | 0.24 | 0.09 | 0.24 | 0.08 | 1.03 | 0.36
CEERS | 0.19 | 0.08 | 0.16 | 0.08 | 1.10 | 0.40
COSMOS | 0.19 | 0.10 | 0.12 | 0.10 | 1.01 | 0.28
UDS | 0.16 | 0.09 | 0.09 | 0.10 | 1.09 | 0.41
SMACS0723 | 0.21 | 0.10 | 0.17 | 0.11 | 1.02 | 0.37
Total | 0.18 | 0.09 | 0.13 | 0.10 | 1.06 | 0.37

Figure 10 displays the distributions of $RFF$ and asymmetry in each redshift
range. We note that there is a positive correlation between $RFF$ and
asymmetry, implying that asymmetric galaxies have high residual fluxes from a
single Sérsic fitting. The $RFF$ distributions consistently show median values
of $RFF\sim 0.13$ at $z<4$, whereas at $z>4$, the median $RFF$ values become
lower than $0.1$. This decrease in $RFF$ at higher redshifts might be due to
the reduced effective radii in both physical and angular sizes of galaxies in
the early universe. In contrast to asymmetry, there are some systematic
variations in $RFF$ values of galaxies depending on the JWST fields, as shown
in  Table 5. These field-to-field variations result from differences in
background fluctuations within the images of the JWST fields. For instance,
the NGDEEP field shows relatively higher $RFF$ with a median value of $0.24$
and $\sigma_{\rm MAD}=0.08$, compared to other fields with a median value of
$0.13$ and $\sigma_{\rm MAD}=0.10$. This is because the JWST images of NGDEEP
have a deeper surface brightness limit, leading to lower background sigma in
Equation 2. However, these systematic variations across JWST fields have a
negligible effect when applying the generous criterion of $RFF>0.5$ to
identify poorly-fit galaxies. In addition, a majority of galaxies classified
as irregulars are selected with the asymmetry criterion of $A>0.32$, so the
$RFF>0.5$ criterion has an insignificant effect on the classification of
irregular galaxies.

Figure 11:  Distributions of Sérsic indices as a function of stellar mass. The
panels are ordered as in  Figure 9. The median values of Sérsic indices in
each stellar mass bin are represented by diamond symbols, along with their
deviations. The black dashed lines mark the criterion of $n=1.5$ to
distinguish disk galaxies and spheroid galaxies. In this figure, we do not
plot median values for stellar mass bins with fewer than five objects.

Figure 11 displays the distribution of Sérsic index as a function of stellar
mass and redshift. Although we applied the $n=1.5$ criterion to classify disk
and spheroid galaxies, our analysis reveal continuous distributions of Sérsic
indices across all redshift ranges. The reason of our choice of $n=1.5$ is
described in Section 3.1 and the Appendix of Park et al. (2022). Overall,
Sérsic indices tend to increase in the high-mass regime of ${\rm
log}~{}M_{\ast}/M_{\odot}\gtrsim 10.5$, indicating that galaxy mass growth via
mergers is accompanied by bulge-like structures with higher Sérsic indices.
Examining the redshift evolution of Sérsic indices, we find that galaxies at
$z>3$ exhibit slightly lower Sérsic index compared to those at $z<3$. With our
classification framework, spheroid galaxies with $n>1.5$ appear more dominant
than disk galaxies in the ${\rm log}~{}M_{\ast}/M_{\odot}\gtrsim 10.5$ range
at $z<3$ and become rare at $z>3$. As shown in  Table 5, the distribution of
Sérsic index is statistically consistent across all the JWST fields with a
very small field-to-field fluctuation of the median $n$.

###  4.2 Morphological Distribution with Stellar Mass and Redshift

Here, we explain the morphological fractions in various stellar mass and
redshift ranges, utilizing data collected from all six JWST fields. It is
worth noting that, unlike the other five blank fields, SMACS0723 has a dense
galaxy cluster at $z=0.39$. This may lead to significant interference from
bright light emitted by cluster members and the intracluster medium. In
addition, the gravitational lensing effect induced by the SMACS0723 cluster
can magnify and distort the morphologies of high-redshift galaxies. Despite
these potential problems, the general trends of morphological fractions from
the small sample of the SMACS0723 field seem to be quite consistent with those
in other fields (see Appendix A). Thus, we use the data from the SMACS0723
field for  Figures 12, 13, and 14, with only excluding galaxies at
$z=0.6-1.2$.

Figure 12:  Morphological fractions as a function of stellar mass in various
redshift bins. The fractions of disks, spheroids, and irregulars are denoted
by blue diamonds, red circles, and green triangles, respectively. We plot the
data from five JWST fields, excluding those from the SMACS0723 field (see
texts). We do not plot the fractions for stellar mass bins with fewer than
five objects.

Figure 12 shows morphological fractions relative to stellar mass in redshift
bins corresponding to different JWST filters. To calculate uncertainties in
each bin, we considered both Poisson noise from the sample size and field-to-
field variations in morphological fractions. Detailed information on
morphological fractions for each JWST field is provided in Appendix A, along
with the calculation method for uncertainties.  Figure 12 reveals a clear
dependence of disk and spheroid fractions on stellar mass. In the ${\rm
log}~{}M_{\ast}/M_{\odot}<10$ mass regime, disk galaxies are dominant
constituting $\gtrsim 60\%$ at all redshift ranges. Spheroid fractions
increase with stellar masses at all redshifts. At $z<3$, spheroid galaxies
become dominant in the high-mass regime, specifically ${\rm
log}~{}M_{\ast}/M_{\odot}>10.5$ at $z<2$ and ${\rm
log}~{}M_{\ast}/M_{\odot}>11$ at $z=2-3$. Beyond $z>3$, spheroid galaxies
consistently show lower proportions than disk galaxies across all stellar mass
ranges. Nevertheless, increasing trends in spheroid fractions still appear
from the mass range of ${\rm log}~{}M_{\ast}/M_{\odot}\gtrsim 10$. In the case
of irregular galaxies, their fractions remain nearly constant at around $20\%$
in all ranges of stellar mass and redshift.

Figure 13:  Morphological fractions as a function of stellar mass (left panel)
and redshift (right panel). Symbols are the same as  Figure 12. Here we also
collect the data from five JWST fields except for the SMACS0723 field.

We plot the combined data of morphological fractions in relation to stellar
mass and redshift in  Figure 13. The left panel shows the dependence of
morphological fraction on stellar mass, combining data at all redshift ranges.
As stellar mass increases, the disk fraction decreases and the spheroid
fraction increases; the crosspoint of disk and spheroid fractions occurs
around ${\rm log}~{}M_{\ast}/M_{\odot}\sim 10.5-11.0$. The irregular fraction
is nearly independent on stellar mass but slightly increases as stellar mass
decreases in the range of ${\rm log}~{}M_{\ast}/M_{\odot}<11$, showing a
larger fraction than spheroid fraction at ${\rm
log}~{}M_{\ast}/M_{\odot}\lesssim 10$. The right panel shows the morphological
fractions with redshifts, providing an insight on cosmic evolution of galaxy
morphology. In the early universe at $z>3$, disk galaxies are a little more
prevalent than at $z<3$, with a fraction of $\gtrsim 70\%$ compared to $\sim
60\%$ at $z<3$, while spheroid fraction increases at this cosmic period. This
is due to the mass dependence of galaxy morphology and the fraction of massive
galaxies increases at lower redshifts. The spheroid fraction is only $\sim
10\%$ at $z>3$, but reaches $\gtrsim 20\%$ at $z<2$. Irregular fraction does
not show a monotonic relation with redshifts. The irregular fraction seems to
reach a maximum during the cosmic noon, $z=2.0-3.0$.

The dominance of disk galaxies at high redshifts and at relatively low stellar
mass implies that the initial morphology of galaxies is disk-like when they
first form. Spheroidal and irregular galaxies also exist beyond $z>3$, but
these galaxies constitute only about $10-20\%$. Changes in morphological
fractions with stellar mass and redshift suggest that spheroid and irregular
galaxies can be formed from initial disk galaxies via mergers or interactions.
In the redshift range of $z\lesssim 2$, some of irregular galaxies seem to
evolve into disk or spheroid galaxies, leading to a decrease in irregular
fraction.

##  5 Discussion

###  5.1 Comparison with the HR5 Results

Figure 14:  Morphological fractions as a function of stellar mass at cosmic
afternoon ($z<1.5$; left panel), cosmic noon ($z=1.5-4.0$; middle panel), and
cosmic morning ($z>4$; right panel). Symbols are the same as  Figure 12. As in
Figure 12, we do not plot the fractions for stellar mass bins with fewer than
five objects.

Our results on the morphology fraction of high-redshift galaxies agree closely
with those of Park et al. (2022) in many ways. It should be noted that Park et
al. (2022)’s prediction had been presented before any JWST results were
announced. In the HR5 simulation, about 70% of the first galaxies with
$9.3<{\rm log}~{}M_{\ast}/M_{\odot}\lesssim 10$ are found to have disk
morphology at $z>4$ (refer to their Figure 11). The rests are divided into
irregulars and spheroids more or less equally. Spheroid and irregular
fractions are nearly independent of redshift and stellar mass within the mass
range, each type constituting around $\sim 1/6$ of the first galaxies. For
comparison with these HR5 results, we depict the morphological fractions with
stellar mass in  Figure 14, categorizing the cosmic epochs of ‘cosmic morning
($z>4$)’, ‘cosmic noon ($z=1.5-4$)’, and ‘cosmic afternoon ($z<1.5$)’,
following Park et al. (2022)’s definition.

The left panel of  Figure 14, showing the morphological fraction of the three
galaxy types in the cosmic morning ($z>4$), illustrates a stunning agreement
between our JWST observational results and HR5 cosmological simulation
predictions. This panel shows that disk galaxies account for about 75% of
galaxies with ${\rm log}~{}M_{\ast}/M_{\odot}\lesssim 10$ at $z>4$, while
spheroid and irregular galaxies account for about $10\%$ and $15\%$,
respectively. These morphological distributions closely resemble those from
the HR5 study, indicating that the theoretical predictions on the morphology
of the first galaxies by HR5 are clearly supported by our JWST results.

In one of the HR5 snapshots during the cosmic noon period ($z\sim 2$), disk
galaxies remain dominant, comprising around $\sim 70\%$ of galaxies in the
${\rm log}~{}M_{\ast}/M_{\odot}\lesssim 10$ range (refer to Figure 20 in Park
et al. (2022)). However, the disk fraction shows a sharp decline in the mass
range ${\rm log}~{}M_{\ast}/M_{\odot}\gtrsim 10.5$. Meanwhile, the spheroid
fraction surpasses disk fraction in the ${\rm log}~{}M_{\ast}/M_{\odot}>11$
range, and irregular fraction also reaches up to $\sim 35\%$ at ${\rm
log}~{}M_{\ast}/M_{\odot}\gtrsim 10.5$. These trends of disk and spheroid
fractions in the HR5 simulation also appear in a very similar fashion in our
JWST analysis, as demonstrated in the middle panel of  Figure 14. The results
from this study show disk dominance up to ${\rm log}~{}M_{\ast}/M_{\odot}\sim
11$ and spheroid dominance at higher masses. However, the irregular fraction
in our observational study shows a somewhat different behavior from that in
the HR5 simulation. The fraction indeed reaches a maximum during the cosmic
noon period, but is nearly constant with $\sim 20\%$ regardless of stellar
mass, in contrast to the jump at ${\rm log}~{}M_{\ast}/M_{\odot}=10-10.5$
found by Park et al. (2022). This discrepancy could arise from differences in
irregular classifiers, such as asymmetry and $RFF$, which may not be matched
perfectly with those used in the HR5 study.

Morphological fractions in the cosmic afternoon ($z<1.5$) in the right panel
of  Figure 14 provide an insight on the morphological evolution of galaxies
from high redshifts ($z>4$) to $z\lesssim 1$. It seems that the disk fraction
in the cosmic afternoon has not changed much since the cosmic noon. On the
other hand, the spheroid fraction shows a noticeable increase, reaching $\sim
70\%$ in the high-mass regime (${\rm log}~{}M_{\ast}/M_{\odot}>11$). The
irregular fraction in the cosmic afternoon exhibits a gentle decline with
increasing stellar mass, in contrast to the constant fraction of $\sim 20\%$
in the cosmic noon. These trends in spheroid and irregular galaxies suggest
that massive irregular galaxies in the cosmic noon or morning are likely to
evolve into spheroid galaxies as they kinematically settle down and maintain
the spheroidal morphology.

In summary, we note the following points from the perspective of these three
cosmic periods. First, disk galaxies are dominant in the cosmic morning with
the fraction of $\sim 75\%$ for those with ${\rm
log}~{}M_{\ast}/M_{\odot}\lesssim 10$, implying that disks are the initial
morphology of the first galaxies. Second, morphological transformations from
disks to spheroids or irregulars occur across all the cosmic epochs,
particularly in the mass range of ${\rm log}~{}M_{\ast}/M_{\odot}\gtrsim
10.5$. Third, spheroid galaxies exist in the cosmic morning at relatively low
masses, but dominates the galaxy population in the cosmic noon and afternoon
in the high-mass regime of ${\rm log}~{}M_{\ast}/M_{\odot}\gtrsim 10.5-11$.
Fourth, the irregular fraction increases from the cosmic morning to the cosmic
noon but decreases again in the cosmic afternoon due to the reduced number of
massive irregular galaxies. Lastly, these trends of morphological fractions
can imply an evolutionary sequence from disks to spheroids or irregulars, with
massive irregular galaxies eventually transforming into spheroids in the
cosmic afternoon. Overall, these findings align very well with the
morphological evolution of galaxies suggested in the HR5 simulation.

###  5.2 Formation and Evolution of Disks from the Early Universe

In this study, we observed a prevalence of disk-type morphology for the
galaxies with stellar masses ${\rm log}~{}M_{\ast}/M_{\odot}\lesssim 10$ at
$z>4$, which aligns with the HR5 results. Interestingly, our findings from
JWST images, capturing the distribution of stellar light, exhibit a good
agreement with the HR5 simulation based on stellar mass distribution. Park et
al. (2022) explained that the emergence of disk-type morphologies in the early
universe can be attributed to the initial angular momentum in protogalactic
clouds. The initial angular momentum is gained from the inflow of cold gas
into protogalaxies along primordial large-scale structures, which is
consistent with the basis of the tidal torque theory (Peebles, 1969; White,
1984). According to this theory, the tidal field on protogalactic regions and
corresponding velocity field are governed by the large-scale dark matter
distribution, which keep galaxies to acquire angular momentum set up by the
initial conditions. The competition between the tidal torque driving galaxies
into disk type and mergers driving galaxies into irregulars and spheroids
determines the morphology distribution. The dominance of disks in the cosmic
morning can be elucidated with these underlying physical processes.

In  Figure 14, our analysis reveals $\sim 20\%$ decrease in disk fraction at
small masses during the cosmic noon and afternoon compared to the cosmic
morning, in contrast to a notable increase in the spheroid fraction in the
high-mass regime. In the low-mass regime (${\rm
log}~{}M_{\ast}/M_{\odot}<10$), there is a slight increase of $5-10\%$ in the
irregular fraction from the cosmic morning to the cosmic noon, indicating a
tendency for disk galaxies in this mass range to transform into irregular
galaxies.

However, it is important to note that mergers do not inevitably lead to
irreversible morphological changes from disks to other types. Despite
undergoing morphological transformations due to mergers or destructive
interactions, the initial disk structures can be restored through the
acquisition of angular momentum via gas accretion from the vorticity-rich
surrounding matter (Welker et al., 2014; Park et al., 2022). Indeed, the HR5
simulation supported this recovery process, demonstrating morphological
fluctuations in initial disk galaxies with changing Sérsic indices (refer to
Figure 12 in Park et al., 2022). Consequently, this disk recovery mechanism
can effectively explain why the disk fraction remains high in the cosmic
morning in spite of high merger rate at corresponding redshifts (Duncan et
al., 2019).

###  5.3 Comparison with Other Observational Studies on Galaxy Morphology

For comparison with other observational studies, it is important to note that
the morphology classification schemes employed in other studies are not
consistent with ours. We applied the parametric classification described in
Section 3.4 for consistent comparison with the HR5 study, but previous studies
have used different morphology classification approaches with various
scientific objectives. Commonly employed methods are visual classification
with measurements of CAS parameters or machine learning techniques applied to
large datasets. Thus, the direct comparison with other studies is difficult.
Furthermore, our parametric classification has a potential problem of
misclassifying the galaxies with parameters around $n\sim 1.5$ or $A\sim 0.32$
as our classification adopts sharp boundaries. Despite these systematic
differences and limitations, we discuss both agreements and disagreements on
morphological distribution and evolution with the results from previous
studies.

Most studies using HST data disagree with our findings of disk dominance in
the early universe. Previous HST studies have suggested a prevalence of
irregular galaxies with disturbed structures at $z\gtrsim 4$, rather than the
disk dominance (Papovich et al., 2005; Cameron et al., 2011; Mortlock et al.,
2013; Huertas-Company et al., 2016). Otherwise, a few studies have
demonstrated high spheroid fractions of $\sim 30\%$ at $z\gtrsim$ in the rest-
frame UV wavelength (Lotz et al., 2006; Dahlen et al., 2007). These
discrepancies with our results may arise from the limited observational
performance of HST for high-redshift galaxy morphology and differences of
observing wavelength ranges in the rest frame. This is also pointed out by
several JWST studies, which suggested that results based on JWST images
exhibit a higher incidence of regular morphologies, such as disks and
spheroids, compared to HST (Ferreira et al., 2022; Jacobs et al., 2023).
Furthermore, JWST images are useful for uncovering hidden disk galaxies within
the noise of HST images (Nelson et al., 2023; Robertson et al., 2023).

Recently, several JWST studies agree with our findings regarding disk
dominance in the early universe. Ferreira et al. (2022, 2023), utilizing
visual classification with JWST fields of the SMACS0723 and CEERS, reported
that disk galaxies occupy $40-60\%$ of the entire galaxy population at
$z=1.5-6.5$. They found that the disk fraction increases by a factor of 10
when using JWST images compared to investigations based solely on HST images.
The spheroid and irregular fractions are nearly constant at $10-20\%$,
$30-40\%$, respectively, with respect to redshifts.

Kartaltepe et al. (2023), also employing visual classification in the CEERS
field, presented a dominant disk fraction of $\sim 60\%$ at $z\sim 3$,
declining to $\sim 30\%$ beyond $z>6$. Similar to our results, Kartaltepe et
al. (2023) noted an increase in spheroid fraction in the high-mass regime
(${\rm log}~{}M_{\ast}/M_{\odot}\gtrsim 10.5$) across all redshift ranges.
They reported pure spheroid and irregular fractions of $\sim 20\%$ at $z>3$,
with a substantial number of galaxies exhibiting mixed morphological
structures due to their classification scheme with seven morphological types.

Huertas-Company et al. (2023) and Tohill et al. (2023) also utilized more
detailed subdivision in their classification scheme than ours, employing CNN
and unsupervised machine learning techniques. They found that irregular and
clumpy galaxies become dominant beyond $z>3$, with high fractions of disk-like
morphologies in ${\rm log}~{}M_{\ast}/M_{\odot}>11$. In summary, previous JWST
studies cannot directly compared with one another as the classification
criteria are all different. We adopt a simple but objective classification
scheme that has been proven to be effective and can be directed compared to
the HR5 cosmological simulation results.

##  6 Summary

In this study, we examine the rest-frame optical morphologies of high-redshift
galaxies using the JWST/NIRCam images obtained from six JWST fields: NEP-TDF,
NGDEEP, CEERS, COSMOS, UDS, and SMACS0723. We select $\sim 18,000$ high-
redshift galaxies with stellar masses of ${\rm log}~{}M_{\ast}/M_{\odot}>9$
and redshifts of $z=0.6-8.0$, derived from the SED fitting procedure. We apply
a parametric morphological classification scheme utilized in the HR5
simulation (Park et al., 2022) and compare our findings with the HR5 results.
The key parameters for morphology classification (Sérsic index, asymmetry, and
$RFF$) were derived from the JWST/NIRCam images through SExtractor photometry
and a single Sérsic fitting with GALFIT. Following the methodology of the HR5
study, our classification scheme categorizes three morphological types: disks
($n<1.5$, $A<0.32$, and $RFF<0.5$), spheroids ($n>1.5$, $A<0.32$, and
$RFF<0.5$), and irregulars ($A>0.32$ or $RFF>0.5$). From these analyses, our
main results can be summarized as follows.

  1. 1.

Our photometric redshift measurements from the SED fitting show an outlier
fraction of $\sim 5\%$ with $\sigma_{\rm NMAD}=0.008$ for objects with
$\chi^{2}_{\rm SED}<20$ and $N_{\rm filt}\geq 6$. This accuracy is superior to
the previous studies using HST-only data, highlighting that JWST photometric
data is beneficial to analyze the properties of a large sample of high-
redshift galaxies. However, several outliers at $z_{\rm spec}<1$ exhibit
overestimated photometric redshifts with $z_{\rm phot}>1$, potentially
indicating low-redshift contaminants in our sample.

  2. 2.

The distribution of asymmetry is strongly concentrated at $A=0.2$ but skewed
toward higher asymmetry values across all redshift ranges. This skewness
occurs due to the existence of asymmetric galaxies with irregular
substructures. The $RFF$ values have positive correlations with asymmetry,
with a median $RFF=0.13$. The Sérsic indices tend to increase with stellar
mass, implying that the formation of bulge-like structures is accompanied with
mass growth of galaxies. The proportion of galaxies with bulge-like structures
($n>1.5$) declines at $z>3$ compared to the $z<3$ universe. Potential biases
of these parameters depending on the JWST fields have negligible effects on
morphology classification.

  3. 3.

Disk galaxies are dominant in the mass range of ${\rm
log}~{}M_{\ast}/M_{\odot}<10.5$, with the disk fraction reaching up to $\sim
80\%$ at $z>4$ in this mass range. However, the disk fraction decreases with
increasing stellar mass across all redshift ranges. This implies that disks
are likely to be the initial morphology of galaxies in the early universe and
experience morphological transformations to spheroids or irregulars.

  4. 4.

Spheroid galaxies are present at $z>6$ with a fraction of $\lesssim 10\%$ of
the galaxy population. However, the spheroid fraction increases as stellar
mass increases and redshift decreases, becoming the dominant type in ${\rm
log}~{}M_{\ast}/M_{\odot}\gtrsim 10.5-11$ at $z<3$. This spheroid dominance in
the high-mass regime implies that the mass growth is associated with
morphological transformation from disks to spheroids.

  5. 5.

Irregular galaxies maintain a relatively constant fraction of $\sim 20\%$ with
respect to stellar mass. In the morphological fraction as a function of
redshift, however, the irregular fraction exhibits a peak at $z\sim 2-3$ and
decreases again below $20\%$ in the $z<1$ universe. This decrease in irregular
fraction is mainly driven by a decrease in the number of massive irregulars
(${\rm log}~{}M_{\ast}/M_{\odot}\gtrsim 10.5$), implying that irregular
galaxies eventually transform into spheroid galaxies in the cosmic afternoon
($z<1.5$).

  6. 6.

The morphological distributions in the cosmic morning ($z>4$) are consistent
with the HR5 simulation, confirming the disk dominance in the early universe.
As suggested in Park et al. (2022), disk galaxies in the cosmic morning evolve
into spheroids and irregulars through merger and gas accretion. Massive
irregular galaxies appear to be an intermediate phase in the morphological
transition from disks to spheroids. These evolutionary tracks of galaxy
morphology align well with the HR5 study, thus demonstrating that this study
effectively verify the HR5 results with the JWST observational data.

The close agreement between our observational morphology distribution
measurement and HR5 simulation results opens the possibility that the origin
and evolution of galaxy morphology can be understood from cosmological
simulations. It remains to be seen how future observational results with
larger number of galaxies and improved redshift accuracy will compare with our
measurement and also with the predictions of future cosmological simulations.

JHL was supported by the National Research Foundation of Korea (NRF) grant
funded by the Korean government (MSIT) (Nos. 2022R1A4A3031306 and
2023R1A2C1006261). HSH acknowledges the support by the National Research
Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No.
2021R1A2C1094577). This work is based on observations made with the
NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the
Mikulski Archive for Space Telescopes at the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST.

##  Appendix A Field-to-field Variations in Morphological Distributions

Figure 15:  Same as the left panel of  Figure 13, but for each JWST field.
Figure 16:  Same as the right panel of  Figure 13, but for each JWST field.

The use of six JWST fields with different depths can introduce variations in
morphological distributions. As shown in  Table 3, the outlier fraction of our
photometric redshift measurements exhibits fluctuations across the JWST
fields, indicating potential variances in the accuracy of stellar masses
derived from the SED fitting. In addition, variations in background
fluctuations can lead to systematic differences in the estimation of key
parameters ($n$, $A$, and $RFF$) for morphology classification. To account for
these variations, we measured the field-to-field deviations of morphological
fractions in the mass and redshift bins, by calculating standard deviations of
the fractions in each bin weighted by the number of samples in each field. The
total uncertainties in  Figures 12, 13, and 14 were determined by summing the
squared field-to-field deviations and Poisson noises of the sample sizes in
the bins. Here we illustrate the morphological fractions in the six JWST
fields as functions of stellar mass and redshift in  Figure 15 and  Figure 16,
respectively.

Figure 15 shows that disk galaxies have a dominant proportion of $\sim 60\%$
in ${\rm log}~{}M_{\ast}/M_{\odot}<10.5$ for all JWST fields. Although disk
fractions generally decrease with increasing stellar mass, this trend is less
discernible in the fields of NEP-TDF and NGDEEP due to their limited sample
sizes in the mass bin of ${\rm log}~{}M_{\ast}/M_{\odot}=11.4-12.0$. Spheroid
fractions increase with stellar mass, but the data in the NGDEEP field suffer
from the small number statistics in the highest mass bin. Irregular fractions
remain constant at $\sim 20\%$, but show large fluctuations in the highest
mass bin due to significant uncertainties. Irregular fractions in the NGDEEP
field are higher than in other fields, due to relatively larger $RFF$ values
resulting from small background sigma values in Equation 2 (see Section 4.1).
Although we exclude the data from the SMACS0723 field due to the concerns
about gravitational lensing effect, the general tendency in this field is
consistent with the stacked data in the left panel of  Figure 13.

In  Figure 16, disk galaxies also maintain dominance across all redshift
ranges in all JWST fields. Decreases in disk fraction from $\sim 70\%$ at
$z\gtrsim 3$ to $\sim 50\%$ at $z\lesssim 3$ are consistently observed in all
fields. Despite fluctuations in the disk fraction at the highest redshift bin
($z=6.0-8.0$), the general trend appears unaffected by field-to-field
variations. Spheroid fractions consistently rise with decreasing redshift,
except for the cases with only five spheroids in the high-redshift bins
($z=4.0-8.0$) in the NGDEEP field. Irregular fractions also exhibit similar
patterns across all JWST fields, peaking at over $\sim 25\%$ around $z\sim
2-3$. Like other morphological types, irregular fractions in the $z=6.0-8.0$
bin are less reliable in the NEP-TDF and NGDEEP fields due to their small
sample sizes (only 3 irregulars in the NEP-TDF and none in the NGDEEP). The
NGDEEP field shows notably higher irregular fractions of about $\sim 40\%$ at
$z=2-3$ compared to other fields, for the reasons mentioned above.

These findings demonstrate that the general trends for each morphological type
are not significantly affected by the field-to-field variations. However,
these variations are still non-negligible, so we reflect these fluctuations
across the JWST fields as the uncertainties of morphological fractions in the
mass and redshift bins.

##  Appendix B Conversion from Three-dimensional to Two-dimensional Asymmetry

Figure 17:  Comparison between the three-dimensional asymmetry ($A_{\rm 3D}$;
x-axis) and the two-dimensional asymmetry ($A_{\rm 2D}$; y-axis). Blue circles
denote the objects with $A_{\rm 3D}>0.4$, which are defined as asymmetric
galaxies in Park et al. (2022), while red circles represent the rests. Note
that the division threshold of $A_{\rm 3D}=0.4$ has been found based on the
actual distribution $A$ of the galaxies identified in HR5. Gray dashed line
marks the asymmetry criterion of $A_{\rm 2D}=0.32$ used in this study.

The asymmetry criterion of $A=0.4$ applied in Park et al. (2022) to
distinguish symmetric galaxies and asymmetric galaxies was established based
on the three-dimensional distribution of stellar particles. Thus, the direct
application of this criterion to the two-dimensional data of JWST images is
not desireable. To address this, we conducted a simple test to determine a
two-dimensional asymmetry criterion that corresponds to the three-dimensional
asymmetry threshold of $A=0.4$ employed based on the actual distribution $A$
of the galaxies identified in HR5.

In  Figure B1, we present a comparison between the three-dimensional asymmetry
($A_{\rm 3D}$) and the two-dimensional asymmetry ($A_{\rm 2D}$) of $2,669$
simulated galaxies at $z\sim 6$ in the HR5 simulation. The two-dimensional
asymmetry was measured using the same methodology described in Section 3.1 of
Park et al. (2022), with the face-on projected stellar mass density profile.
We searched for the optimal threshold value of $A_{\rm 2D}>X$ corresponding to
$A_{\rm 3D}>0.4$, which minimizes contamination and loss of completeness
simultaneously.

Contamination quantifies the presence of contaminated sources within
asymmetric objects selected with the $A_{\rm 2D}>X$ criterion, which was
calculated as the number fraction of objects with $A_{\rm 2D}>X$ and $A_{\rm
3D}<0.4$ relative to those with $A_{\rm 2D}>X$ and $A_{\rm 3D}>0.4$.
Completeness represents the fraction of objects with $A_{\rm 2D}>X$ and
$A>0.4$ relative to the total number of those with $A_{\rm 3D}>0.4$, so loss
of completeness is the opposite probability of the completeness. We determined
the threshold for two-dimensional asymmetry at $X=0.32$ (denoted by the black
dashed line in the figure), yielding minimum values for both contamination and
loss of completeness about $\sim 10\%$. This two-dimensional asymmetry
criterion of $A_{\rm 2D}=0.32$ was subsequently applied in our morphology
classification, as detailed in Section 3.4.

## References

  * Abraham et al. (1996) Abraham, R. G., van den Bergh, S., Glazebrook, K., et al. 1996, ApJS, 107, 1. doi:10.1086/192352
  * Abraham & van den Bergh (2001) Abraham, R. G. & van den Bergh, S. 2001, Science, 293, 1273. doi:10.1126/science.1060855
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33. doi:10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123. doi:10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167. doi:10.3847/1538-4357/ac7c74
  * Bagley et al. (2023) Bagley, M. B., Pirzkal, N., Finkelstein, S. L., et al. 2023, arXiv:2302.05466. doi:10.48550/arXiv.2302.05466
  * Bamford et al. (2009) Bamford, S. P., Nichol, R. C., Baldry, I. K., et al. 2009, MNRAS, 393, 1324. doi:10.1111/j.1365-2966.2008.14252.x
  * Barden et al. (2012) Barden, M., Häußler, B., Peng, C. Y., et al. 2012, MNRAS, 422, 449. doi:10.1111/j.1365-2966.2012.20619.x
  * Benítez (2000) Benítez, N. 2000, ApJ, 536, 571. doi:10.1086/308947
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393. doi:10.1051/aas:1996164
  * Bezanson et al. (2022) Bezanson, R., Labbe, I., Whitaker, K. E., et al. 2022, arXiv:2212.04026. doi:10.48550/arXiv.2212.04026
  * Bradley et al. (2023) Bradley, L. D., Coe, D., Brammer, G., et al. 2023, ApJ, 955, 13. doi:10.3847/1538-4357/acecfe
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503. doi:10.1086/591786
  * Brammer (2023a) Brammer, Gabriel. (2023a). grizli (1.8). Zenodo. doi:10.5281/zenodo.7701739
  * Brammer (2023b) Brammer, Gabriel. (2023b). eazy-py (0.6.7). Zenodo. doi:10.5281/zenodo.8268031
  * Buitrago et al. (2013) Buitrago, F., Trujillo, I., Conselice, C. J., et al. 2013, MNRAS, 428, 1460. doi:10.1093/mnras/sts124
  * Calvi et al. (2012) Calvi, R., Poggianti, B. M., Fasano, G., et al. 2012, MNRAS, 419, L14. doi:10.1111/j.1745-3933.2011.01168.x10.1086/122140
  * Cameron et al. (2011) Cameron, E., Carollo, C. M., Oesch, P. A., et al. 2011, ApJ, 743, 146. doi:10.1088/0004-637X/743/2/146
  * Carnall et al. (2023) Carnall, A. C., Begley, R., McLeod, D. J., et al. 2023, MNRAS, 518, L45. doi:10.1093/mnrasl/slac136
  * Casey et al. (2023) Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31. doi:10.3847/1538-4357/acc2bc
  * Coe et al. (2019) Coe, D., Salmon, B., Bradač, M., et al. 2019, ApJ, 884, 85. doi:10.3847/1538-4357/ab412b
  * Conselice et al. (2000) Conselice, C. J., Bershady, M. A., & Jangren, A. 2000, ApJ, 529, 886. doi:10.1086/308300
  * Conselice (2003) Conselice, C. J. 2003, ApJS, 147, 1. doi:10.1086/375001
  * Conselice et al. (2005) Conselice, C. J., Blackburne, J. A., & Papovich, C. 2005, ApJ, 620, 564. doi:10.1086/426102
  * Conselice et al. (2008) Conselice, C. J., Rajgor, S., & Myers, R. 2008, MNRAS, 386, 909. doi:10.1111/j.1365-2966.2008.13069.x
  * Dahlen et al. (2007) Dahlen, T., Mobasher, B., Dickinson, M., et al. 2007, ApJ, 654, 172. doi:10.1086/508854
  * Dressler (1980) Dressler, A. 1980, ApJ, 236, 351. doi:10.1086/157753
  * Duncan et al. (2019) Duncan, K., Conselice, C. J., Mundy, C., et al. 2019, ApJ, 876, 110. doi:10.3847/1538-4357/ab148a
  * Dunlop et al. (2021) Dunlop, J. S., Abraham, R. G., Ashby, M. L. N., et al. 2021, JWST Proposal. Cycle 1, 1837
  * Fasano et al. (2015) Fasano, G., Poggianti, B. M., Bettoni, D., et al. 2015, MNRAS, 449, 3927. doi:10.1093/mnras/stv500
  * Ferreira et al. (2022) Ferreira, L., Adams, N., Conselice, C. J., et al. 2022, ApJ, 938, L2. doi:10.3847/2041-8213/ac947c
  * Ferreira et al. (2023) Ferreira, L., Conselice, C. J., Sazonova, E., et al. 2023, ApJ, 955, 94. doi:10.3847/1538-4357/acec76
  * Finkelstein et al. (2017) Finkelstein, S. L., Dickinson, M., Ferguson, H. C., et al. 2017, JWST Proposal ID 1345. Cycle 0 Early Release Science, 1345
  * Finkelstein et al. (2021) Finkelstein, S. L., Papovich, C., Pirzkal, N., et al. 2021, JWST Proposal. Cycle 1, 2079
  * Finkelstein et al. (2023) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13. doi:10.3847/2041-8213/acade4
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1. doi:10.1051/0004-6361/202039657
  * Galametz et al. (2013) Galametz, A., Grazian, A., Fontana, A., et al. 2013, ApJS, 206, 10. doi:10.1088/0067-0049/206/2/10
  * Green (2018) Green, G. M. 2018, The Journal of Open Source Software, 3, 695. doi:10.21105/joss.00695
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357. doi:10.1038/s41586-020-2649-2
  * Hoyos et al. (2011) Hoyos, C., den Brok, M., Verdoes Kleijn, G., et al. 2011, MNRAS, 411, 2439. doi:10.1111/j.1365-2966.2010.17855.x
  * Hoyos et al. (2012) Hoyos, C., Aragón-Salamanca, A., Gray, M. E., et al. 2012, MNRAS, 419, 2703. doi:10.1111/j.1365-2966.2011.19918.x
  * Hubble (1926) Hubble, E. P. 1926, ApJ, 64, 321. doi:10.1086/143018
  * Huertas-Company et al. (2009) Huertas-Company, M., Tasca, L., Rouan, D., et al. 2009, A&A, 497, 743. doi:10.1051/0004-6361/200811255
  * Huertas-Company et al. (2016) Huertas-Company, M., Bernardi, M., Pérez-González, P. G., et al. 2016, MNRAS, 462, 4495. doi:10.1093/mnras/stw1866
  * Huertas-Company et al. (2023) Huertas-Company, M., Iyer, K. G., Angeloudi, E., et al. 2023, arXiv:2305.02478. doi:10.48550/arXiv.2305.02478
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90. doi:10.1109/MCSE.2007.55
  * Hwang & Park (2009) Hwang, H. S. & Park, C. 2009, ApJ, 700, 791. doi:10.1088/0004-637X/700/1/791
  * Ikeda et al. (2023) Ikeda, R., Morishita, T., Tsukui, T., et al. 2023, MNRAS, 523, 6310. doi:10.1093/mnras/stad1692
  * Jacobs et al. (2023) Jacobs, C., Glazebrook, K., Calabrò, A., et al. 2023, ApJ, 948, L13. doi:10.3847/2041-8213/accd6d
  * Kartaltepe et al. (2021) Kartaltepe, J., Casey, C. M., Bagley, M., et al. 2021, JWST Proposal. Cycle 1, 1727
  * Kartaltepe et al. (2023) Kartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2023, ApJ, 946, L15. doi:10.3847/2041-8213/acad01
  * Kron (1980) Kron, R. G. 1980, ApJS, 43, 305. doi:10.1086/190669
  * Lee et al. (2021) Lee, J., Shin, J., Snaith, O. N., et al. 2021, ApJ, 908, 11. doi:10.3847/1538-4357/abd08b
  * Lotz et al. (2004) Lotz, J. M., Primack, J., & Madau, P. 2004, AJ, 128, 163. doi:10.1086/421849
  * Lotz et al. (2006) Lotz, J. M., Madau, P., Giavalisco, M., et al. 2006, ApJ, 636, 592. doi:10.1086/497950
  * Margalef-Bentabol et al. (2016) Margalef-Bentabol, B., Conselice, C. J., Mortlock, A., et al. 2016, MNRAS, 461, 2728. doi:10.1093/mnras/stw1451
  * Mortlock et al. (2013) Mortlock, A., Conselice, C. J., Hartley, W. G., et al. 2013, MNRAS, 433, 1185. doi:10.1093/mnras/stt793
  * Nelson et al. (2023) Nelson, E. J., Suess, K. A., Bezanson, R., et al. 2023, ApJ, 948, L18. doi:10.3847/2041-8213/acc1e1
  * Noirot et al. (2023) Noirot, G., Desprez, G., Asada, Y., et al. 2023, MNRAS, 525, 1867. doi:10.1093/mnras/stad1019
  * Ormerod et al. (2023) Ormerod, K., Conselice, C. J., Adams, N. J., et al. 2023, arXiv:2309.04377. doi:10.48550/arXiv.2309.04377
  * Papovich et al. (2005) Papovich, C., Dickinson, M., Giavalisco, M., et al. 2005, ApJ, 631, 101. doi:10.1086/429120
  * Park & Choi (2005) Park, C. & Choi, Y.-Y. 2005, ApJ, 635, L29. doi:10.1086/499243
  * Park et al. (2007) Park, C., Choi, Y.-Y., Vogeley, M. S., Gott, J. R., & Blanton, M. R. 2007, ApJ, 658, 898.
  * Park et al. (2008) Park, C., Gott, J. R., & Choi, Y.-Y. 2008, ApJ, 674, 784. doi:10.1086/524192
  * Park & Choi (2009) Park, C. & Choi, Y.-Y. 2009, ApJ, 691, 1828. doi:10.1088/0004-637X/691/2/1828
  * Park & Hwang (2009) Park, C. & Hwang, H. S. 2009, ApJ, 699, 1595. doi:10.1088/0004-637X/699/2/1595
  * Park et al. (2022) Park, C., Lee, J., Kim, J., et al. 2022, ApJ, 937, 15. doi:10.3847/1538-4357/ac85b5
  * Peebles (1969) Peebles, P. J. E. 1969, ApJ, 155, 393. doi:10.1086/149876
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., et al. 2002, AJ, 124, 266. doi:10.1086/340952
  * Peng et al. (2010) Peng, C. Y., Ho, L. C., Impey, C. D., et al. 2010, AJ, 139, 2097. doi:10.1088/0004-6256/139/6/2097
  * Perrin et al. (2014) Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, Proc. SPIE, 9143, 91433X. doi:10.1117/12.2056689
  * Petrosian (1976) Petrosian, V. 1976, ApJ, 210, L53. doi:10.1086/182301
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13. doi:10.1051/0004-6361/201525830
  * Pontoppidan et al. (2022) Pontoppidan, K. M., Barrientes, J., Blome, C., et al. 2022, ApJ, 936, L14. doi:10.3847/2041-8213/ac8a4e
  * Postman et al. (2005) Postman, M., Franx, M., Cross, N. J. G., et al. 2005, ApJ, 623, 721. doi:10.1086/428881
  * Rigby et al. (2023) Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001. doi:10.1088/1538-3873/acb293
  * Robertson et al. (2023) Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, ApJ, 942, L42. doi:10.3847/2041-8213/aca086
  * Sandage (1961) Sandage, A. 1961, Washington: Carnegie Institution, 1961
  * Schlawin et al. (2020) Schlawin, E., Leisenring, J., Misselt, K., et al. 2020, AJ, 160, 231. doi:10.3847/1538-3881/abb811
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525. doi:10.1086/305772
  * Tohill et al. (2023) Tohill, C.-B., Bamford, S., Conselice, C., et al. 2023, arXiv:2306.17225. doi:10.48550/arXiv.2306.17225
  * Valentino et al. (2023) Valentino, F., Brammer, G., Gould, K. M. L., et al. 2023, ApJ, 947, 20. doi:10.3847/1538-4357/acbefa
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261. doi:10.1038/s41592-019-0686-2
  * Vulcani et al. (2011) Vulcani, B., Poggianti, B. M., Aragón-Salamanca, A., et al. 2011, MNRAS, 412, 246. doi:10.1111/j.1365-2966.2010.17904.x
  * Ward et al. (2023) Ward, E. M., de la Vega, A., Mobasher, B., et al. 2023, arXiv:2311.02162. doi:10.48550/arXiv.2311.02162
  * Weaver et al. (2023) Weaver, J. R., Cutler, S. E., Pan, R., et al. 2023, arXiv:2301.02671. doi:10.48550/arXiv.2301.02671
  * Welker et al. (2014) Welker, C., Devriendt, J., Dubois, Y., et al. 2014, MNRAS, 445, L46. doi:10.1093/mnrasl/slu106
  * White (1984) White, S. D. M. 1984, ApJ, 286, 38. doi:10.1086/162573
  * Windhorst et al. (2022) Windhorst, R. A., Ashcraft, T. A., Broadhurst, T. J., et al. 2022, JWST Proposal. Cycle 1, 2738
  * Windhorst et al. (2023) Windhorst, R. A., Cohen, S. H., Jansen, R. A., et al. 2023, AJ, 165, 13. doi:10.3847/1538-3881/aca163