# Spatially Resolved Stellar Populations of $0.3<z<6.0$ Galaxies in
WHL0137$-$08 and MACS0647+70 Clusters as Revealed by JWST: How do Galaxies
Grow and Quench Over Cosmic Time?

Abdurro’uf Center for Astrophysical Sciences, Department of Physics and
Astronomy, The Johns Hopkins University, 3400 N Charles St. Baltimore, MD
21218, USA  Space Telescope Science Institute (STScI), 3700 San Martin Drive,
Baltimore, MD 21218, USA  Dan Coe Space Telescope Science Institute (STScI),
3700 San Martin Drive, Baltimore, MD 21218, USA  Association of Universities
for Research in Astronomy (AURA), Inc. for the European Space Agency (ESA)
Center for Astrophysical Sciences, Department of Physics and Astronomy, The
Johns Hopkins University, 3400 N Charles St. Baltimore, MD 21218, USA  Intae
Jung Space Telescope Science Institute (STScI), 3700 San Martin Drive,
Baltimore, MD 21218, USA  Henry C. Ferguson Space Telescope Science Institute
(STScI), 3700 San Martin Drive, Baltimore, MD 21218, USA  Gabriel Brammer
Cosmic Dawn Center (DAWN), Copenhagen, Denmark  Niels Bohr Institute,
University of Copenhagen, Jagtvej 128, Copenhagen, Denmark  Kartheik G. Iyer
Hubble Fellow  Columbia Astrophysics Laboratory, Columbia University, 550 West
120th Street, New York, NY 10027, USA  Larry D. Bradley Space Telescope
Science Institute (STScI), 3700 San Martin Drive, Baltimore, MD 21218, USA
Pratika Dayal Kapteyn Astronomical Institute, University of Groningen, P.O.
Box 800, 9700 AV Groningen, The Netherlands  Rogier A. Windhorst School of
Earth and Space Exploration, Arizona State University, Tempe, AZ 85287-1404,
USA  Adi Zitrin Physics Department, Ben-Gurion University of the Negev, P.O.
Box 653, Be’er-Sheva 84105, Israel  Ashish Kumar Meena Physics Department,
Ben-Gurion University of the Negev, P.O. Box 653, Be’er-Sheva 84105, Israel
Masamune Oguri Center for Frontier Science, Chiba University, 1-33 Yayoi-cho,
Inage-ku, Chiba 263-8522, Japan  Department of Physics, Graduate School of
Science, Chiba University, 1-33 Yayoi-Cho, Inage-Ku, Chiba 263-8522, Japan
Jose M. Diego Instituto de Física de Cantabria (CSIC-UC). Avda. Los Castros
s/n. 39005 Santander, Spain  Vasily Kokorev Kapteyn Astronomical Institute,
University of Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands
Paola Dimauro INAF - Osservatorio Astronomico di Roma, via di Frascati 33,
00078 Monte Porzio Catone, Italy  Angela Adamo Department of Astronomy, Oskar
Klein Centre, Stockholm University, AlbaNova University Centre, SE-106 91
Stockholm, Sweden  Christopher J. Conselice Jodrell Bank Centre for
Astrophysics, University of Manchester, Oxford Road, Manchester UK  Brian
Welch Department of Astronomy, University of Maryland, College Park, MD 20742,
USA  Observational Cosmology Lab, NASA Goddard Space Flight Center, Greenbelt,
MD 20771, USA  Center for Research and Exploration in Space Science and
Technology, NASA/GSFC, Greenbelt, MD 20771  Eros Vanzella INAF – OAS,
Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Gobetti
93/3, I-40129 Bologna, Italy  Tiger Yu-Yang Hsiao Center for Astrophysical
Sciences, Department of Physics and Astronomy, The Johns Hopkins University,
3400 N Charles St. Baltimore, MD 21218, USA  Xinfeng Xu Center for
Astrophysical Sciences, Department of Physics and Astronomy, The Johns Hopkins
University, 3400 N Charles St. Baltimore, MD 21218, USA  Namrata Roy Center
for Astrophysical Sciences, Department of Physics and Astronomy, The Johns
Hopkins University, 3400 N Charles St. Baltimore, MD 21218, USA  Celia R.
Mulcahey Center for Astrophysical Sciences, Department of Physics and
Astronomy, The Johns Hopkins University, 3400 N Charles St. Baltimore, MD
21218, USA

###### Abstract

We study the spatially resolved stellar populations of 444 galaxies at
$0.3<z<6.0$ in two clusters (WHL0137-08 and MACS0647+70) and a blank field,
combining imaging data from HST and JWST to perform spatially resolved
spectral energy distribution (SED) modeling using piXedfit. The high spatial
resolution of the imaging data combined with magnification from gravitational
lensing in the cluster fields allows us to resolve a large fraction of our
galaxies (109) to sub-kpc scales. At redshifts around cosmic noon and higher
($2.5\lesssim z\lesssim 6.0$), we find mass doubling times to be independent
of radius, inferred from flat specific star formation rate (sSFR) radial
profiles and similarities between the half-mass and half-SFR radii. At lower
redshifts ($1.5\lesssim z\lesssim 2.5$), a significant fraction of our star-
forming galaxies show evidence for nuclear starbursts, inferred from centrally
elevated sSFR, and a much smaller half-SFR radius compared to the half-mass
radius. At later epochs, we find more galaxies suppress star formation in
their centers but are still actively forming stars in the disk. Overall, these
trends point toward a picture of inside-out galaxy growth consistent with
theoretical models and simulations. We also observe a tight relationship
between the central mass surface density and global stellar mass with $\sim
0.38$ dex scatter. Our analysis demonstrates the potential of spatially
resolved SED analysis with JWST data. Future analysis with larger samples will
be able to further explore the assembly of galaxy mass and the growth of their
structures.

Galaxy evolution (594) — Galaxy formation (595) – Galaxy clusters (584) –
Galaxy quenching (2040)

††journal: ApJ††facilities: HST, JWST††software: astropy (Astropy
Collaboration et al., 2013, 2018, 2022), piXedfit (Abdurro’uf et al., 2021,
2022c), SExtractor (Bertin & Arnouts, 1996), sep (Barbary, 2016), photutils
(Bradley et al., 2022a), grizli (Brammer et al., 2022), eazypy (Brammer et
al., 2008).

##  1 Introduction

Over the last few decades, multiwavelength studies of galaxies throughout
cosmic history reveal that the global star formation rate density (SFRD) in
the universe was increasing with cosmic time from the reionization epoch and
reached a peak at $z\sim 2$ ($\sim 3.5$ Gyr after the Big Bang; cosmic noon)
after which it declined exponentially toward the present day (Madau &
Dickinson, 2014). In this picture, it is estimated that $\sim$25% of the
present-day stellar mass density (SMD) was formed before the peak of the
cosmic SFRD, around half of the SMD was formed during $0.7<z<2.0$, and another
$\sim$25% was formed since $z=0.7$ (i.e., around the last half of the
universe’s age; Madau & Dickinson, 2014). Although the cosmic SFRD at early
cosmic time is still debated due to the dust obscuration effects (see e.g.,
Fudamoto et al. 2021; Casey et al. 2018), an emerging picture is that cosmic
SMD increases with cosmic time since the epoch of reionization, which is
believed to take place before $z\sim 6$ (e.g., Treu et al., 2013; McGreer et
al., 2015; Dayal & Ferrara, 2018).

Observations also revealed that most of the star formation occurs in galaxies
that lie in the so-called star-forming main sequence (SFMS), which is a tight
nearly linear correlation between the integrated (i.e., global) star formation
rate (SFR) and stellar mass ($M_{*}$; Brinchmann et al., 2004; Daddi et al.,
2007; Elbaz et al., 2007; Noeske et al., 2007; Whitaker et al., 2012, 2014;
Speagle et al., 2014; Salmon et al., 2015; Tomczak et al., 2016; Santini et
al., 2017; Iyer et al., 2018; Leslie et al., 2020; Leja et al., 2022; Popesso
et al., 2023). This relation has been shown to hold at any epoch with a nearly
constant scatter ($\sim 0.3$ dex; Whitaker et al., 2012; Speagle et al., 2014;
Popesso et al., 2023), suggesting that galaxies grow in mass over cosmic time
in a state of self-regulated semi-equilibrium (e.g., Bouché et al., 2010;
Daddi et al., 2010; Genzel et al., 2010; Tacchella et al., 2016b, 2020).
Understanding this process in detail and the mechanisms that shut down star
formation in galaxies and move them out of the SFMS onto the “quenched”
population requires knowledge of not only integrated galaxy properties but
also spatially resolved structures within galaxies.

The study of spatially resolved properties of galaxies with integral field
spectroscopy (IFS) and high spatial resolution imaging has been providing
important insights toward a better understanding of galaxy evolution. Among
the important findings is the realization that some of the well-known scaling
relations observed on global scales originated from similar relations on kpc
scales within galaxies (see a review by Sánchez, 2020). This includes the
spatially resolved star-forming main sequence (rSFMS), a relationship between
the SFR surface density ($\Sigma_{\rm SFR}$) and $M_{*}$ surface density
($\Sigma_{*}$), a local analog of the global SFMS (e.g., Sánchez et al., 2013;
Wuyts et al., 2013; Cano-Díaz et al., 2016; Abdurro’uf & Akiyama, 2017; Hsieh
et al., 2017; Abdurro’uf & Akiyama, 2018; Lin et al., 2019; Enia et al.,
2020). Recent studies found that the rSFMS relation (and hence the global
SFMS) originated from two more fundamental relations on kpc scales: the
resolved Kennicutt–Schmidt ($\Sigma_{\rm SFR}$ versus $\rm{H}_{2}$ mass
surface density; $\Sigma_{\rm{H}_{2}}$) and molecular gas main-sequence
relations ($\Sigma_{*}$ versus $\Sigma_{\rm{H}_{2}}$) (e.g., Lin et al., 2019;
Baker et al., 2022; Abdurro’uf et al., 2022b; Morselli et al., 2020). This
emphasizes the necessity of studying the spatially resolved properties of
galaxies.

Spatially resolved studies of high redshift galaxies ($z\sim 1-4$) have hinted
at how galaxies assembled their structures. The emerging picture from these
studies is that galaxies grow their mass in an inside-to-outside manner (i.e.,
inside-out growth scenario; e.g., van Dokkum et al., 2013; Nelson et al.,
2012; Morishita et al., 2015; Nelson et al., 2016) and cease their star
formation activities in a similar manner (i.e., inside-out quenching scenario;
e.g., Tacchella et al., 2015; Jung et al., 2017; Abdurro’uf & Akiyama, 2018;
Tacchella et al., 2018; Ellison et al., 2018; Bluck et al., 2020). Nelson et
al. (2016) analyzed the spatially resolved distributions (on kpc scales) of
$\text{H}\alpha$ emission and stellar mass of $0.7<z<1.5$ galaxies using the
Hubble Space Telescope (HST)/WFC3 grism data from the 3D-HST survey (Skelton
et al., 2014). They found that the spatial distribution of $H_{\alpha}$
emission in the galaxies is more extended than the stellar mass distribution,
suggesting that the past star formation in the galaxies has accumulated
stellar mass in the center and now the star formation progresses outward to
assemble the disk. Tacchella et al. (2015) analyzed the spatial distributions
of SFR and $M_{*}$ of $\sim 30$ star-forming galaxies at $z\sim 2$ using IFS
data from the SINS/zC-SINF survey (Förster Schreiber et al., 2018). They
observed that massive galaxies ($M_{*}\gtrsim 10^{11}M_{\odot}$) in their
sample have a centrally-suppressed specific SFR (sSFR) radial profile and a
massive central spheroid that is as dense as the centers of local early-type
galaxies. In contrast to this, less massive galaxies in their sample have
broadly flat sSFR radial profiles. This trend indicates that massive galaxies
at this epoch might have started a quenching process in their central regions
and assembled a mature bulge.

The buildup of the central stellar mass density is likely correlated with the
quenching process in galaxies. The central stellar mass density within a 1 kpc
radius ($\Sigma_{*,1\rm{kpc}}$) has been shown to be a good predictor for
quiescence, where galaxies with high $\Sigma_{*,1\rm{kpc}}$ tend to be red and
quiescent, whereas galaxies with low $\Sigma_{*,1\rm{kpc}}$ tend to be blue
and star-forming (e.g., Fang et al., 2013; Tacchella et al., 2015, 2016a;
Barro et al., 2017; Jung et al., 2017; Whitaker et al., 2017; Dimauro et al.,
2022). It has also been shown that $\Sigma_{*,1\rm{kpc}}$ is tightly
correlated with the global $M_{*}$, suggesting that $M_{*}$ of galaxies grow
hand-in-hand with the central mass density. In this
$\Sigma_{*,1\rm{kpc}}$–$M_{*}$ relation, quiescent galaxies reside in a
sequence at the tip of the overall relationship and have a shallower slope
than the relation with star-forming galaxies only, indicating a formation of a
matured bulge in the quiescent galaxies (Fang et al., 2013; Tacchella et al.,
2015; Barro et al., 2017).

Galaxies also grow their sizes hand-in-hand with the global $M_{*}$, as
indicated by the size–mass relation (e.g., Shen et al., 2003; van der Wel et
al., 2014; Suess et al., 2019). Previous studies have shown that star-forming
and quiescent galaxies follow very different size–mass relations where
quiescent galaxies tend to be more compact (i.e., having smaller size) in all
$M_{*}$ and exhibit steeper relation than the star-forming galaxies (van der
Wel et al., 2014; Yang et al., 2021). A possible explanation for this trend is
that star-forming galaxies build their mass at all radii by mostly in-situ
star formation, whereas quiescent galaxies grow inside-out through mergers
(e.g., van Dokkum et al., 2015).

Previous studies, some of which are mentioned above, have used HST for
resolving galaxies out to $z\sim 3$, roughly a limit where galaxies can be
resolved well by the telescope, given its spatial resolution and depth.
Furthermore, the wavelength coverage of HST only covers the rest-frame
ultraviolet (UV) and a small portion of the optical at $z\sim 3$, making it
difficult to robustly derive $M_{*}$ as well as the other stellar population
properties, which typically requires a rest-frame near-infrared (NIR). Forcing
to include NIR imaging from the ground-based telescopes would need to
sacrifice the spatial resolution of HST (e.g., Jung et al., 2017). With the
advent of the James Webb Space Telescope (JWST) NIRCam observations (Rieke et
al., 2022; Rigby et al., 2022), with its high spatial resolution, depth, and
its coverage in NIR, now we can push the analysis of spatially resolved SED of
galaxies to higher redshifts. Some very recent studies have used JWST/NIRCam
imaging to study the internal structures and morphology of galaxies at $z>3$
(e.g., Ferreira et al., 2022; Chen et al., 2022; Kartaltepe et al., 2022;
Pérez-González et al., 2022; Giménez-Arteaga et al., 2022), and even resolving
a lensed galaxy at $z\sim 11$ (Hsiao et al., 2022).

In this paper, we use imaging data from HST/ACS and JWST/NIRCam to analyze the
spatially resolved SEDs of $0.3<z<6.0$ galaxies in the sightlines of
WHLJ013719.8–082841 (hereafter WHL0137$-$08; RA = 01:37:25.0, DEC =
$-$08:27:23, J2000; $z=0.566$; Wen et al. 2012; Wen & Han 2015) and
MACSJ0647.7+7015 (hereafter MACS0647+70; RA = 06:47:50.03, DEC =
$+$70:14:49.7, J2000; $z=0.591$; Ebeling et al. 2007) clusters and examine the
spatial distributions of their stellar populations. Our main goal is to get
hints on the assembly of galaxy structures over cosmic time, especially how
galaxies build their stellar masses and quench their star formation
activities. The high spatial resolution of JWST/NIRCam combined with
magnification from gravitational lensing in the cluster fields, allows us to
resolve high-redshift galaxies down to sub-kpc scales. Our method using
piXedfit (Abdurro’uf et al., 2022c) can simultaneously process imaging data,
perform pixel binning to optimize the signal-to-noise (S/N) ratio of the
spatially resolved SEDs, and perform SED fitting. The wavelength coverage of
HST/ACS and JWST/NIRCam allow us to get full coverage of the rest-frame UV to
NIR for the majority of our sample, which can give a strong constraint on
model SEDs and break the age–dust–metallicity degeneracy (see Appendix B).
While IFS observation at $z\gtrsim 2$ is lacking, our analysis in this paper
provides a good alternative for the analysis of spatially resolved SED of high
redshift galaxies. Our analysis in this paper is one of the first robust
spatially resolved SED analyses of hundreds of galaxies using JWST data.
Abdurro’uf et al. (2021) have demonstrated the capabilities of spatially
resolved SED fitting using piXedfit on local galaxies. In particular, it gives
robust SFR on kpc scales when rest-frame UV–NIR photometry is available, which
is consistent with the SFR derived from $\text{H}\alpha$ emission maps (dust-
corrected based on the Balmer decrement) from the MaNGA IFS survey (Bundy et
al., 2015).

The paper is organized as follows. In Section 2, we present the data and
sample galaxies. We describe the spatially resolved SED fitting methodology in
Section 3 and present our results in Section 4, which include the radial
profiles of some key stellar population properties, comparison between the
compactness of the spatial distributions of SFR and $M_{*}$, and
$\Sigma_{*,1\rm{kpc}}$–$M_{*}$ relation. In Section 5, we further discuss our
results, focusing on the evolutionary trends with redshift and the
implications to the study of galaxy evolution.

Throughout this paper, we assume the Chabrier (2003) initial mass function
(IMF) with a mass range of $0.1-100M_{\odot}$ and cosmological parameters of
$\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and
$H_{0}=70\text{km}\text{s}^{-1}\text{Mpc}^{-1}$.

##  2 Data and Sample

###  2.1 Observational Data

Table 1: HST and JWST Imaging Data Used in the Spatially Resolved SED Fitting

Telescope | Camera | Filter | Wavelength | Deptha | PSF FWHMb
---|---|---|---|---|---
|  |  |  | WHL0137$-$08 | MACS0647+70 | WHL0137$-$08 | MACS0647+70
|  |  | (µm) | (AB mag) | (AB mag) | (arcsec) | (arcsec)
HST | ACS/WFC | F435W | 0.37–0.47 | 27.7 | 28.0 | 0.11 | 0.11
HST | ACS/WFC | F475W | 0.4–0.55 | 28.5 | 28.2 | 0.11 | 0.11
HST | ACS/WFC | F555W | 0.46–0.62 |  $\cdots$ | 28.7 |  $\cdots$ | 0.11
HST | ACS/WFC | F606W | 0.47–0.7 | 28.3 | 28.3 | 0.11 | 0.11
HST | ACS/WFC | F625W | 0.54–0.71 |  $\cdots$ | 27.9 |  $\cdots$ | 0.11
HST | ACS/WFC | F775W | 0.68–0.86 |  $\cdots$ | 27.8 |  $\cdots$ | 0.08
HST | ACS/WFC | F814W | 0.7–0.95 | 28.7 | 28.5 | 0.11 | 0.11
JWST | NIRCam | F090W | 0.8–1.0 | 28.3 |  $\cdots$ | 0.04 |  $\cdots$
JWST | NIRCam | F115W | 1.0–1.3 | 28.4 | 28.1 | 0.04 | 0.04
JWST | NIRCam | F150W | 1.3–1.7 | 28.5 | 28.3 | 0.06 | 0.06
JWST | NIRCam | F200W | 1.7–2.2 | 28.7 | 28.4 | 0.06 | 0.06
JWST | NIRCam | F277W | 2.4–3.1 | 29.1 | 28.9 | 0.11 | 0.11
JWST | NIRCam | F356W | 3.1–4.0 | 29.3 | 29.0 | 0.11 | 0.11
JWST | NIRCam | F410M | 3.8–4.3 | 28.6 |  $\cdots$ | 0.16 |  $\cdots$
JWST | NIRCam | F444W | 3.8–5.0 | 29.0 | 28.8 | 0.16 | 0.16

Note. — a$5\sigma$ point source AB magnitude limit measured within a $0\farcs
2$ diameter circular aperture. bPSF FWHM is based on empirical measurement as
described in Appendix C.

####  2.1.1 JWST Observations

We obtain JWST/NIRCam imaging data of WHL0137$-$08 cluster from Cycle 1
General Observers (GO) 2282 program (PI Coe) and MACS0647+70 cluster from GO
1433 program (PI Coe). The WHL0137$-$08 cluster was observed in July 2022,
while the MACS0647+70 cluster was observed on 23 September 2022. The GO 2282
program aims at further investigating Earendel (Welch et al., 2022a, b) and
the Sunrise Arc (Vanzella et al., 2022). The JWST/NIRCam data from this
program consist of eight filters (F090W, F115W, F150W, F200W, F277W, F365W,
F410M, and F444W) spanning a wavelength range of $0.8-5.0~{}\micron$. The GO
1433 program is intended to observe the triply-lensed galaxy MACS0647–JD at
$z\sim 11$ (Coe et al., 2013; Hsiao et al., 2022). This program obtained
JWST/NIRCam imaging in six filters (F115W, F150W, F200W, F277W, F365W, and
F444W) spanning $1-5~{}\micron$. The exposure time of each filter in the two
programs is $2104$ seconds. It achieves $5\sigma$ limiting AB magnitude of
$28.0$ to $29.0$ in a $r=0\farcs 2$ diameter circular aperture.

For each filter, we obtained four dithers using INTRAMODULEBOX primary dithers
to cover the $4-5\arcsec$ gap between the sort wavelength (SW;
$\lambda<2.4\micron$) detectors, improve the spatial resolution of final
drizzled images, and minimize the impact of image artifacts and bad pixels. In
each observation, we obtained NIRCam imaging over two $2\farcm 26\times
2\farcm 26$ fields separated by 40$\farcm$5, covering a total area of 10.2
arcmin2. In the observation of WHL0137$-$08 cluster, the NIRCam module B was
centered at the cluster while module A covered a nearby field centered $\sim$
2$\farcm$9 from the cluster center (hereafter called “blank field”). On the
other hand, the MACS0647+70 cluster was centered at module A, and module B
observes a blank field nearby to it.

####  2.1.2 HST Data

We obtain HST imaging data of the WHL0137$-$08 cluster from the Reionization
Lensing Cluster Survey (RELICS) HST Treasury program (GO 14096; Coe et al.,
2019). The RELICS program obtained the first HST imaging of the WHL0137$-$08
cluster in 2016 with three orbits of ACS (F435W, F606W, and F814W) and two
orbits of WFC3/IR (F105W, F125W, F140W, and F160W) data spanning
$0.4-1.7~{}\micron$. Two follow-up HST imaging programs (GO 15842 and GO
16668; PI: Coe) have thus far obtained an additional 5 orbits of HST ACS
imaging in F814W, 2 orbits in F475W, and 4 orbits with WFC3/IR in F110W.

The HST imaging data of the MACS0647+70 cluster are taken from multiple
programs. Overall, MACS0647+70 has been observed in total of 39 orbits of HST
imaging in 17 filters. The cluster was first observed by programs GO 9722 (PI
Ebeling) and GO 10493, 10793 (PI Gal-Yam) in the ACS F555W and F814W filters.
Then additional imaging in 15 filters (WFC3/UVIS, ACS, and WFC3/IR, spanning
$0.2-1.7~{}\micron$ ) was obtained by the Cluster Lensing and Supernova Survey
with Hubble (CLASH; Postman et al. 2012; GO 12101, PI Postman). Finally,
additional imaging in WFC3/IR F140W was obtained as part of a grism
spectroscopy program (GO 13317, PI Coe).

It is important to note that the nearby blank fields to the WHL0137$-$08 and
MACS0647+70 clusters that were observed with NIRCam are not covered in the HST
observations described above. In this work, we analyze galaxies in three
fields: WHL0137$-$08 cluster field, MACS0647+70 cluster, and the NIRCam blank
field nearby to the WHL0137$-$08 (hereafter simply called blank field). We do
not analyze galaxies in the NIRCam blank field of MACS0647+70 because it is
observed in fewer filters than the blank field of WHL0137$-$08 and it does not
have F090W observation, which prevents us from selecting galaxies at $z<2$ in
this field as their photometry do not cover the rest-frame $4000$ Å break. For
the WHL0137$-$08 and the blank field, we use 4 HST/ACS filters (F435W, F475W,
F606W, and F814W) and 8 JWST/NIRCam filters, whereas for the MACS0647+70, we
use 7 HST/ACS filters (F435W, F475W, F555W, F606W, F625W, F775W, and F814W)
and 6 JWST/NIRCam filters. We do not use HST/WFC3 IR filters to get high
spatial resolution possible while still getting sufficiently wide wavelength
coverage with the HST/ACS and JWST/NIRCam. Please refer to Table 1 for
information on limiting magnitudes and the point spread function (PSF) sizes
of our data.

###  2.2 Sample Galaxies

We use grizli v4 photometric catalogs (will be described in Section 3.1) to
select our sample galaxies in the three fields (WHL0137$-$08, blank field, and
MACS0647+70). The catalogs provide aperture fluxes and photometric redshifts
with which we select our sample. The sample selection is described in the
following. First, we select galaxies that have integrated signal-to-noise
(S/N) ratio $>5$ in all JWST filters that are available for the fields. This
is to ensure that we will have galaxies with good photometry in at least JWST
filters. This initial cut selects 1322 (out of 2718), 1278 (out of 3032), and
1331 (out of 2660) galaxies in the WHL0137$-$08, blank field, and MACS0647+70,
respectively. We do not apply the same S/N criteria on HST filters because it
would exclude more galaxies as they have lower S/N than JWST filters.

We further cut the sample galaxies based on their redshift to get a sufficient
coverage of the rest-frame UV–NIR. For galaxies in the WHL0137$-$08 and
MACS0647+70, which are observed by both JWST and HST, we select galaxies at
$0.3<z<6.0$, whereas for galaxies in the blank field, which do not have HST
observations, we select galaxies at $1.3<z<6.0$. This redshift cut ensures
that the rest-frame $4000$ Å break is covered. This cut further reduces the
sample to 1258, 581, and 1257 galaxies in the WHL0137$-$08, blank field, and
MACS0647+70, respectively. After that, we do a visual inspection to exclude
galaxies that appear to be very small (i.e., unresolved) and in a merger
(i.e., one segmentation region having multiple cores or multiple galaxies in
one segmentation region, despite possible interlopers). This further reduces
the sample to 354, 239, and 220 galaxies in the WHL0137$-$08, blank field, and
MACS0647+70, respectively.

We perform spatially resolved SED analysis on the galaxies in this initial
sample. A detailed description of the methodology will be given in Section 3.
Once the analysis is done, we inspect the results of all the galaxies and
further exclude galaxies that seem to have bad SED fitting results based on
the $\chi^{2}$ values of the fitting to the integrated SEDs within the central
effective radius (see Section 3.5) and the average $\chi^{2}$ values of the
fitting to the first 20 spatial bins (see Section 3.4 for the definition of
the spatial bin). We exclude galaxies that have $\chi^{2}>20$ for the SEDs
within the central effective radius and average $\chi^{2}>40$ for the first 20
spatial bins. We note that $\chi^{2}$ value can be unrealistically high if
systematic uncertainties of the photometry are not properly accounted for.
Besides this, there is still uncertainty around the zero-point calibration of
NIRCam photometry in the current early observations (e.g., Boyer et al., 2022;
Finkelstein et al., 2022). Therefore, we visually inspect the SED fitting
results of each galaxy using similar plots as shown in Figure 4. We find that
in most cases, NIRCam fluxes are well fitted by our models, better than
HST/ACS fluxes. This might be due to the shallower depths (and lower S/N) of
HST compared to JWST. The $\chi^{2}$ values above are high enough to get a
sufficient number of galaxies and low enough to get good quality SED fitting
results. This results in our final sample, consisting of 243, 91, and 110
galaxies in the WHL0137$-$08, blank field, and MACS0647+70, respectively.
Figure 1 shows the distributions of redshifts and $M_{*}$ of our sample
galaxies.

We note that our sample selection may possibly bias toward selecting
relatively massive, bright, and resolved (i.e., big) galaxies in each
redshift. However, due to the lensing magnification in the cluster fields, we
expect to detect on average lower-mass galaxies with better spatial resolution
than in the normal fields. The small number of galaxies and the limited volume
sampled might make our sample to be not representative of the general
population of galaxies. However, since we do not make inferences on the
average trends or number densities as a function of global properties (e.g.,
$M_{*}$), but instead we show trends in individual galaxies, our results still
provide useful insights on the evolution of galaxy structures. We also ignore
the possible contamination by the Active Galactic Nucleus (AGN) host galaxies
in our current study because of the lack of diagnostics for identifying them
using our current data.

Figure 1: Distributions of $M_{*}$ and redshifts of the sample galaxies
analyzed in this work, which consists of galaxies in two cluster fields
(WHL0137$-$08 and MACS0647+70) and a blank field (a nearby field centered
$\sim$ 2$\farcm$9 from the WHL0137$-$08 that was observed with JWST/NIRCam).
The global $M_{*}$ shown here are derived by summing up $M_{*}$ of pixels (in
the galaxy’s region) obtained from our spatially resolved SED fitting.

##  3 Methodology

###  3.1 Data Reduction and Photometric Catalog

We use the grizli pipeline (Brammer et al., 2022) to process the HST FLT and
the JWST pipeline-calibrated level-2 imaging data. The JWST data were
processed using the calibration pipeline v1.5.3 with CRDS context
jwst_0942.pmap, which includes photometric calibrations based on in-flight
data. The JWST level-2 imaging data were then scaled with detector-dependent
factors (Brammer, 2022) based on a NIRCam flux calibration using the standard
star J1743045. Our photometric zeropoints described here are similar to those
obtained by the JWST Resolved Stellar Populations ERS program (Boyer et al.,
2022; Nardiello et al., 2022) who analyzed the M92 globular cluster. We also
check the consistency of our calibration with the more recent one based on CAL
program data jwst_0989.pmap and find out that they are consistent within 3% in
all filters analyzed here.

In processing the JWST data, the grizli pipeline applies a correction to
reduce the effect of $1/f$ noise and masks “snowballs”111https://jwst-
docs.stsci.edu/data-artifacts-and-features/snowballs-artifact effect caused by
the large cosmic ray impacts to the NIRCam detectors. Besides this, the grizli
pipeline also corrects for the “wisps”222https://jwst-docs.stsci.edu/jwst-
near-infrared-camera/nircam-features-and-caveats/nircam-claws-and-wisps, which
is a faint, diffuse stray light features that appear at the same detector
locations in NIRCam images and most prominent in the A3, B3, and B4 detectors
in the F150W and F200W images.

The grizli pipeline aligns the HST and JWST imaging data to a common world
coordinate system which is registered based on the GAIA DR3 catalogs (Gaia
Collaboration et al., 2021). The images are then drizzled to a common pixel
grid using the astrodrizzle (Koekemoer et al., 2003; Hoffmann et al., 2021).
The 17 HST filters and 4 JWST NIRCam long-wavelength (LW) filters (F277W,
F356W, F410M, and F444W) are drizzled to a spatial sampling of 0$\farcs$04 per
pixel while the JWST short-wavelength (SW) filters (F090W, F115W, F150W, and
F200W) are drizzled to a spatial sampling of 0$\farcs$02 per pixel.

Source detection is then performed on a weighted sum of the drizzled NIRCam
images in all filters using sep (Barbary, 2016; Bertin & Arnouts, 1996).
Fluxes are then calculated for each source in three circular apertures,
$0\farcs 36$, $0\farcs 5$, and $0\farcs 7$. Then photometric redshift
measurement is performed using the $0\farcs 5$ aperture SEDs employing
eazypy(Brammer et al., 2008). This code fits observed photometry using a set
of templates added in a non-negative linear combination. The processed imaging
data along with the photometric catalog are publicly
available333https://cosmic-spring.github.io/index.html. These data products
have also been used in some recent studies (Welch et al., 2022b; Bradley et
al., 2022b; Hsiao et al., 2022; Vanzella et al., 2022; Meena et al., 2022).

###  3.2 Analysis of Post-processed Imaging Data

In this work, we combine the post-processed HST and JWST imaging data (in up
to 13 filters) into a common spatial resolution (i.e., PSF size) and sampling
(i.e., pixel size) for extracting the spatially resolved SEDs of our sample
galaxies. These spatially resolved SEDs are then fitted with models to infer
the underlying properties of the stellar populations. We use
piXedfit444https://github.com/aabdurrouf/piXedfit (Abdurro’uf et al., 2021,
2022c) throughout this analysis. Basically, this process includes three main
tasks: image processing, pixel binning, and SED fitting. We will briefly
describe these steps in the following.

The image processing is carried out automatically using piXedfit. For each
galaxy, we first crop stamp images with a size of 6$\farcs$04 $\times$
6$\farcs$04 (corresponding to 302$\times$302 pixels in NIRCam SW and
151$\times$151 pixels in NIRCam LW and HST/ACS filters) centered at the
galaxy. We then perform background subtraction to each stamp image using
photutils (Bradley et al., 2022a). Next, we perform point spread function
(PSF) matching to homogenize the spatial resolution across filters. We degrade
the spatial resolution of the images to match the resolution of the F444W
filter, which has the lowest spatial resolution (see Table 1). For this, we
generate the empirical PSFs of HST/ACS and JWST/NIRcam filters along with the
convolution kernels using photutils package (see Appendix C). The PSF matching
is carried out by convolving the stamp images with the convolution kernels.
After PSF matching, we register all the stamp images to a common spatial
sampling of 0$\farcs$04 per pixel. At the end, we have multiband stamp images
with a size of 151$\times$151 pixels for each galaxy in our sample.

###  3.3 Constructing Photometric Data Cubes

piXedfit further processes the stamp images to produce photometric data cubes.
First, it defines a galaxy’s region of interest. For each galaxy, segmentation
maps are first produced in all filters using sep (Barbary, 2016) and those
maps are then merged together into a single map. In the segmentation process,
we use the same parameters for all filters as follows. We set the detection
threshold (thresh), the number of thresholds for deblending (deblend_nthresh),
and the minimum contrast ratio for deblending (deblend_cont) to be $2.0$,
$40$, and $0.005$, respectively.

In some cases, the merged segmentation map is larger than expected, as can be
inferred from the maps of multiband fluxes. This can be caused by some
factors, for example, interference from neighboring objects that are not
separated well by the deblending process. We visually inspect the merged
segmentation map of each galaxy to find out this issue. To deal with this, we
tweak the deblending parameters to get cleaner segmentation maps or ignore the
segmentation map in some filters that have this deblending issue, then merge
them again.

Once the galaxy’s region is defined, then the fluxes of pixels within the
region are calculated. We use the PHOTFLAM keyword in the header of the grizli
imaging data products to convert the pixel value into flux density in the
units of erg $\text{s}^{-1}\text{cm}^{-2}$Å-1. The data cubes are then stored
in FITS files. Figure 2 shows examples of the maps of multiband fluxes of
galaxies in three fields analyzed in this work. The color images shown in the
leftmost panels are created using Trilogy555https://github.com/dancoe/trilogy
(Coe et al., 2012).

Figure 2: Examples of the maps of multiband fluxes produced from the image
processing. An example of one galaxy is shown for each field, from top to
bottom: WHL0137$-$08 (observed in 12 filters), blank field (8 filters), and
MACS0647+70 (13 filters). The galaxy ID is based on our grizli v4 public
catalog.

###  3.4 Pixel binning

The SEDs of pixels are usually noisy and might not provide sufficient
constraint to the models if the fitting is performed to them. Therefore, we
perform pixel binning using piXedfit to optimize the signal-to-noise ratio
(S/N) of the spatially resolved SEDs. Basically, this process bins neighboring
pixels to achieve a certain S/N ratio threshold that can be set in multiple
bands. The unique pixel binning scheme in piXedfit, which takes into account
the similarity in SED shape among pixels, allows for achieving a sufficient
S/N ratio in multiple filters of interest while preserving important spatial
information at the pixel level. A detailed description of this pixel binning
scheme is given in Abdurro’uf et al. (2021).

We assume the following parameters in the pixel binning process. We refer the
reader to Abdurro’uf et al. (2021) for more information about the parameters.
We set S/N thresholds to $5$ in all JWST NIRCam filters. We do not set the S/N
threshold to HST filters because the S/N ratio of pixels in the HST images is
low, especially for galaxies at high redshifts. Setting an S/N threshold on
the HST filters would put a strong constraint in the pixel binning process
which can produce a coarser binning map and loosing important spatial
information from the original images.

The rest of the binning parameters are as follows. We set a minimum diameter
of 7 pixels, which is larger than the PSF FWHM size of our data cubes, a
reduced $\chi^{2}$ limit of 5 in the evaluation of the similarity of SED
shape. We refer to F277W flux in determining the brightest pixel to be the
center of a spatial bin. We store new data cubes produced from this pixel
binning process into FITS files. The total number of spatial bins in our
sample galaxies is 24999. Figure 3 shows examples of the pixel binning maps
produced from this process.

Figure 3: Examples of pixel binning results. The pixel binning process
achieves a minimum S/N ratio of 5 in all JWST NIRCam filters.

###  3.5 Spatially Resolved SED Fitting

Table 2: Free Parameters in the SED Modeling and the Assumed Priors.

Parameter | Description | Prior | Sampling/Scale
---|---|---|---
$M_{*}$ |  Stellar mass |  Uniform: min$=\log(s_{\rm best})-2$, max$=\log(s_{\rm best})+2$a | Logarithmic
$Z_{*}$ |  Stellar metallicity |  Uniform: min$=-2.0+\log(Z_{\odot})$, max$=0.2+\log(Z_{\odot})$ | Logarithmic
$t$ |  Time since the onset of star formation ($\text{age}_{\rm sys}$) b |  Uniform: min$=-1.0$, max = age of the universe at the galaxy’s redshift | Logarithmic
$\tau$ |  Parameter that controls the peak time in the double power-law SFH modelb |  Uniform: min$=-1.5$, max$=1.14$ | Logarithmic
$\alpha$ |  Parameter in the double power-law SFH model that controls the slope of the falling star formation episodeb |  Uniform: min$=-2.0$, max$=2.0$ | Logarithmic
$\beta$ |  Parameter in the double power-law SFH model that controls the slope of the rising star formation episodeb |  Uniform: min$=-2.0$, max$=2.0$ | Logarithmic
$\hat{\tau}_{1}$ |  Dust optical depth of the birth cloud in the Charlot & Fall (2000) dust attenuation law |  Uniform: min$=0.0$, max$=4.0$ | Linear
$\hat{\tau}_{2}$ |  Dust optical depth of the diffuse ISM in the Charlot & Fall (2000) dust attenuation law |  Uniform: min$=0.0$, max$=4.0$ | Linear
$n$ |  Power law index in the Charlot & Fall (2000) dust attenuation law |  Uniform: min$=-2.2$, max$=0.4$ | Linear

Note. — a$s_{\rm best}$ is the normalization of model SED derived from the
initial fitting with the $\chi^{2}$ minimization method (see Section 4.2.1 in
Abdurro’uf et al. 2021). bThe mathematical form of the double power-law SFH is
given in Abdurro’uf et al. (2021, Equation 7 therein).

Once we have the binned data cubes, we perform SED fitting to the SEDs of
individual spatial bins in our sample galaxies. Here we use the SED fitting
module in piXedfit. The SED fitting in piXedfit uses a fully Bayesian
technique. We refer the reader to Abdurro’uf et al. (2021) for a detailed
description of the SED modeling and fitting methods as well as comprehensive
tests of its capabilities. In Appendix A, we perform SED-fitting tests using
mock SEDs to demonstrate the robustness of our SED fitting method on combined
HST and JWST photometry. Moreover, in Appendix B we discuss how NIRCam
photometry can potentially help in breaking the degeneracies among age, dust,
and metallicity in SED fitting. In the following, we provide a brief
description of the method and some assumptions applied in our SED fitting.

We use the Flexible Stellar Population Synthesis code (FSPS; Conroy et al.,
2009; Conroy & Gunn, 2010). It includes the nebular emission modeling that
uses the CLOUDY code (Ferland et al., 1998, 2013). In this work, we assume the
Chabrier (2003) initial mass function (IMF), Padova isochrones (Girardi et
al., 2000; Marigo & Girardi, 2007; Marigo et al., 2008), and MILES stellar
spectral library (Sánchez-Blázquez et al., 2006; Falcón-Barroso et al., 2011).
For the star formation history model, we assume an analytic model in the form
of a double power-law. It has been shown in Abdurro’uf et al. (2021) that this
SFH form can give robust inferences of the stellar population properties and
even SFH of galaxies, as tested using synthetic SEDs of simulated galaxies in
the IllustrisTNG simulations. For simulating the effect of dust attenuation,
we use the two-component dust attenuation law of Charlot & Fall (2000). This
dust attenuation law gives an extra attenuation to stars younger than $10$
Myr, in addition to standard attenuation in the diffuse ISM. We model the
attenuation due to the intergalactic medium using Inoue et al. (2014) model.
Since we do not have photometric data that covers the rest-frames mid-infrared
(MIR) and far-infrared (FIR), we switch off the modeling of dust emission and
AGN dusty torus emission in the analysis throughout this work. The SED
modeling has 9 free parameters. We summarize these parameters along with the
assumed priors in Table 2. We assume a constant ionization parameter ($U$) of
$0.01$ in the modeling of the nebular emission.

In the current analysis, we rely on photometric redshift for all of our sample
galaxies because we do not have spectroscopic observations at the moment we
carry out this analysis. To get redshift estimates of the galaxies, we perform
SED fitting with piXedfit in which redshift is let to be free in the fitting.
For this, we fit integrated SED within the effective radius of the galaxies.
The effective radius is measured in F444W image stamp using GALFIT(Peng et al.
2002; see Section 4.2). This is performed to get SEDs with high S/N while
reducing contamination from noisy SEDs of pixels in the outskirt regions. In
this fitting, we apply a prior on redshift in the form of a Gaussian function
centered at the photometric redshift estimated by the eazypy taken from the
grizli catalog (see Section 3.1). We set a width of $0.5$ for this Gaussian
prior. This fitting is performed to derive redshift only. We then use this
redshift information for the SED fitting of all spatial bins in the galaxy, in
which we fix the redshift. We apply the Markov Chain Monte Carlo (MCMC) method
in piXedfit. In the SED fitting for redshift determination, we set the number
of walkers to 100 and the number of steps per walker to 1000. For the SED
fitting of spatial bins, we use 100 walkers and less number of steps per
walker (600) for reducing computational time.

We show examples of SED fitting results of two galaxies in Figure 4, one
galaxy from the WHL0137$-$08 cluster field (top panel) and the other galaxy
from the blank field (bottom panel). For each galaxy, we show best-fit SEDs in
the top right panel. The observed and best-fit photometric SEDs are shown with
square and circle symbols, respectively. The SED in black represents the
integrated SED within the effective radius, while those in other colors are
for 5 examples of spatial bins in the galaxies. The corner plot in the bottom
left side shows the posterior probability distribution functions (PPDF) of the
model parameters obtained from the fitting on the integrated SED within the
effective radius. Above this corner plot, we show the PDFs of $M_{*}$ and SFR
of the example spatial bins. The best-fit spectra shown in the plot are drawn
from the MCMC sampler chains. Therefore, it is possible to get a slight shift
in wavelength between the best-fit spectra of the central SED (where $z$ is
free in the fitting) and that of the spatial bins (where $z$ is fixed in the
fitting). This wavelength shift reflects the uncertainty of the estimated
redshift. Finally, in the bottom right panel, we show the maps of stellar
population properties, including the $M_{*}$ surface density ($\Sigma_{*}$),
SFR surface density ($\Sigma_{\rm SFR}$), mass-weighted age, $A_{V,1}$
($\equiv 1.086\times\hat{\tau}_{1}$), $A_{V,2}$ ($\equiv
1.086\times\hat{\tau}_{2}$), and metallicity.

Figure 4: Examples of SED fits of a galaxy in the WHL0137$-$08 cluster (top
panel) and the blank field (bottom panel). The SED plots show the best-fit
SEDs from the fitting to integrated SED within the effective radius (black
color) and 5 examples of spatial bins (in colors). The corner plots show the
posterior probability distributions of the parameters obtained from fitting to
the central integrated SED. Above this corner plot, we show PDFs of $M_{*}$
and SFR of the 5 example spatial bins. Finally, the maps of stellar population
properties derived from this analysis are shown in the bottom right corner.

###  3.6 Lens Modeling

To estimate the magnifications due to the gravitational lensing effect by the
clusters, we use the lens models constructed by our team. For the WHL0137$-$08
cluster, we use the same lens models that were used for analyzing the Earendel
and the Sunrise Arc in Welch et al. (2022a), which were made publicly
available666https://relics.stsci.edu/lens_models/outgoing/whl0137-08/. These
lens models were generated using four independent lens modeling software
packages: Light-Traces-Mass (LTM, Zitrin et al., 2009, 2015; Broadhurst et
al., 2005), Glafic (Oguri, 2010), WSLAP (Diego et al., 2005, 2007), and
Lenstool (Kneib et al., 1993; Jullo et al., 2007; Jullo & Kneib, 2009). Please
refer to Welch et al. (2022a) for detailed information about each model.
Sample galaxies located in the blank field are expected to have only weak
magnifications of $\mu\lesssim 1.1$. With the multiple lens models available
for this cluster, we estimate the total and tangential (i.e., linear)
magnifications ($\mu$ and $\mu_{t}$, respectively) of each galaxy by taking
the average values. In this way, we account for the modeling uncertainties.
Based on the standard deviation values, we find that the magnifications do not
vary a lot among the models. The median standard deviations of $\mu$ and
$\mu_{t}$ are $0.11$ and $0.10$ dex, respectively.

The lens models for MACS0647+70 cluster have been constructed in the past
using the HST imaging data. The first lens model for MACS0647+70, before the
CLASH survey, was provided by Zitrin et al. (2011) using LTM method. With the
addition of HST imaging data from CLASH, new lens models were established
using various methods, including Lenstool, LTM, WSLAP, and LensPerfect (Coe et
al., 2008). These lens models have been used in previous studies in CLASH
(e.g., Coe et al., 2013; Zitrin et al., 2015; Chan et al., 2017). Now with the
addition of JWST NIRCam imaging data, which add on many new strongly-lensed
multiple-image candidates (thanks to its high spatial resolution and depth), a
new lens model has been established using the dPIEeNFW method (Zitrin et al.,
2015) with some modifications. Detailed information on this lens modeling of
MACS0647+70 cluster along with the list of the multiple-image systems
considered in the modeling are given in Meena et al. (2022). This new method
has also been implemented in several clusters using JWST NIRCam data (Pascale
et al., 2022; Roberts-Borsani et al., 2022; Hsiao et al., 2022; Williams et
al., 2022). We used this lens model constructed by Meena et al. (2022) for
galaxies in the MACS0647+70 field. We correct the $M_{*}$ and SFR obtained
from SED fitting for the lensing magnification by dividing them with $\mu$. We
also correct size or radius measurements by dividing them with $\mu_{t}$.

##  4 Results

###  4.1 Integrated Properties

Before analyzing the spatially resolved properties of our sample galaxies, we
first present their integrated (i.e., global) properties. To bring it into the
context of the global demographics of galaxies, we plot our sample on the
integrated star-forming main sequence (SFMS) diagram, as shown in the left
panel of Figure 5. The integrated $M_{*}$ and SFR of a galaxy are derived by
summing up the values in pixels obtained from the spatially resolved SED
fitting. Due to our limited sample, we plot all our galaxies on the SFMS
diagram instead of dividing them into a number of redshift bins and examining
the SFMS relation in each bin. This can cause a broad distribution as shown in
the figure. Different symbols represent the fields where the galaxies are
located (WHL0137$-$08, blank field, and MACS0647+70), whereas color-coding
represents redshift grouping, where we divide the redshift range into five
bins. The dashed lines show the SFMS relations at the median redshifts of the
five redshift bins, calculated using the prescription from Speagle et al.
(2014). The lines are colored based on the redshift groups.

We then classify our sample galaxies into star-forming, green valley, and
quiescent groups based on their positions with respect to the SFMS ridge line
at the redshift of the galaxies. We define star-forming, green-valley, and
quiescent galaxies as those having $\text{SFR}>\text{SFR}_{\rm
MS}(z,M_{*})-0.4$ dex, $\text{SFR}_{\rm
MS}(z,M_{*})-0.4\geq\text{SFR}>\text{SFR}_{\rm MS}(z,M_{*})-1.0$ dex, and
$\text{SFR}\leq\text{SFR}_{\rm MS}(z,M_{*})-1.0$ dex, respectively, where
$\text{SFR}_{\rm MS}(z,M_{*})$ is the SFMS ridge line for exact $z$ and
$M_{*}$ of the individual galaxies. With this selection criteria, we have 219,
108, and 117 total numbers of the star-forming, green-valley, and quiescent
galaxies from the three fields, respectively. We will use these classified
samples throughout the analysis in this paper to investigate the differences
in spatially resolved properties of galaxies in various evolutionary stages.
The right panel of Figure 5 show the distributions of these galaxy groups on
the SFMS diagram.

Figure 5: Left panel: Integrated (i.e., global) $M_{*}$ and SFR of our sample
galaxies. Different symbols represent different fields where the galaxies are
located, whereas the color-coding represents redshift grouping. The dashed
lines show the global SFMS relations at the median redshifts of the 5 redshift
groups calculated using the prescription from Speagle et al. (2014). The lines
are colored based on the redshift groups. Right panel: Distribution of the
star-forming, green-valley, and quiescent galaxies in our sample that are
classified as having $\text{SFR}>\text{SFR}_{\rm MS}(z,M_{*})-0.4$ dex,
$\text{SFR}_{\rm MS}(z,M_{*})-0.4\geq\text{SFR}>\text{SFR}_{\rm
MS}(z,M_{*})-1.0$ dex, and $\text{SFR}\leq\text{SFR}_{\rm MS}(z,M_{*})-1.0$
dex, respectively. $\text{SFR}_{\rm MS}(z,M_{*})$ is the SFMS ridge line that
is calculated for exact $z$ and $M_{*}$ of the individual galaxies.

To get a sense of how the global specific SFR (sSFR$\equiv\text{SFR}/M_{*}$)
evolves with cosmic time in our sample galaxies, we plot the sSFR against
redshift in Figure 6. We can see a clear trend of decreasing global sSFR with
cosmic time and an increasing number of quiescent galaxies along the way. In
our sample, quiescent galaxies start to emerge from $z\sim 3$, around $2$ Gyr
after the Big Bang. To compare our global sSFR trend with the similar trend
from previous studies, we plot sSFR$(z)$ inferred from the SFMS normalization
based on Speagle et al. (2014) prescription. We calculate sSFR$(z)$ with four
$M_{*}$ of $10^{8.5}$, $10^{10}$, $10^{11}$, and $10^{12}$ $M_{\odot}$ and
show them in the figure as black dashed lines. We see an overall agreement
between the evolutionary trend of sSFR in our sample and that expected based
on the evolution of the SFMS normalization. We also show the global sSFR of
local ($z\sim 0$) galaxies from Abdurro’uf & Akiyama (2017) and Abdurro’uf et
al. (2022a) (personal communication), which were derived from spatially
resolved SED fitting. Abdurro’uf & Akiyama (2017) analyzed 93 spiral galaxies
at $0.01<z<0.02$ using imaging data from the GALEX (Morrissey et al., 2007)
and SDSS (York et al., 2000) surveys. Abdurro’uf et al. (2022a) applied
piXedfit for analyzing 10 nearby galaxies using imaging data in more than 20
filters spanning from Far-ultraviolet (FUV) to FIR.

Previous studies have classified passive galaxies using various methods. One
of the methods is by comparing the Hubble time ($t_{H}$) with the mass
doubling time (i.e., inverse of sSFR). Basically, this method defines
quiescent galaxies as those having $\text{sSFR}<1/t_{H}$. The red dashed line
in Figure 6 represents $1/t_{H}$. We can see that our quiescent galaxies lie
below this line, indicating that our classification method is consistent with
that based on $t_{H}$.

Figure 6: Evolution of the integrated specific SFR (sSFR) with redshift. The
black dashed lines represent sSFR evolution of SFMS galaxies with
$M_{*}=10^{8.5}$, $10^{10}$, $10^{11}$, and $10^{12}$ $M_{\odot}$ (decreasing
normalization) as inferred from the normalization of the SFMS, calculated
using the prescription from Speagle et al. (2014). The red dashed line
represents $1/t_{H}$ where $t_{H}$ is the Hubble time. Almost all our
quiescent galaxies lie below this line, indicating that their mass doubling
timescale is longer than $t_{H}$.

###  4.2 Radial Profiles of the Stellar Population Properties

As we have shown the global properties of the sample galaxies and classified
them into star-forming, green-valley, and quiescent groups, now we will
analyze their spatially resolved properties. We start by presenting the radial
profiles of the stellar population properties to get a sense of how the
properties vary radially within the galaxies. To derive the radial profiles,
first, we perform 2D single-component Sérsic fitting using GALFIT (Peng et
al., 2002) on F444W stamp image of each galaxy to get their ellipticities,
position angles, and central coordinates. We then use this information to
define elliptical annuli in the radial profile calculation. The radial
profiles are derived from the 2D maps of properties obtained from the
spatially resolved SED fitting by averaging values of pixels within the
annuli. Since galaxies have a wide range of size, we normalize the radius by
the half-mass radius ($R_{e}$), which is the radius that covers half of the
integrated $M_{*}$. We use radial increment ($\delta r$) of $0.3R_{e}$. Thanks
to the gravitational lensing effect, we can resolve many of our galaxies down
to sub-kpc scales (109 galaxies in our sample have delensed $R_{e}<1$ kpc).

Figure 7 shows the radial profiles of the stellar mass surface density
($\Sigma_{*}$). We divide the sample into 5 bins of redshift and 4 bins of
$M_{*}$ to see how the radial profiles vary with global $M_{*}$ and cosmic
time. Moreover, we indicate the star-forming, green valley, and quiescent
galaxies with different colors, in a similar way as in Figure 6. For groups
that contain at least 5 galaxies, we show average radial profiles with tick
lines. Some interesting trends from Figure 7 are the following. At each
redshift bin, more massive galaxies tend to have higher normalization of
$\Sigma_{*}(r)$ than less massive ones, indicating that the excess in mass
happens across the entire radius. Moreover, we also see that quiescent
galaxies tend to have higher $\Sigma_{*}(r)$ normalization than the star-
forming and green-valley galaxies in all redshift. This is especially clear in
the most massive groups. It is also interesting to see that $\Sigma_{*}(r)$
profiles have a negative gradient (i.e., decreasing mass with increasing
radius) in all redshift and mass bins, although the profiles seem to be
shallower at higher redshifts.

Figure 7: Radial profiles of the stellar mass surface density ($\Sigma_{*}$).
The sample galaxies are divided into 4 bins of global $M_{*}$ and 5 bins of
redshift. At each group, we further classify the galaxies into star-forming,
green-valley, and quiescent groups and indicate them with different colors.
For sub-groups that contain at least 5 galaxies, we show average radial
profiles with tick lines. At each redshift bin, more massive galaxies tend to
have higher $\Sigma_{*}(r)$ normalization than less massive galaxies,
indicating that the excess in mass happens across the galaxy region. Quiescent
galaxies tend to have higher $\Sigma_{*}(r)$ normalization than star-forming
in all redshifts. This is especially clear in high $M_{*}$ bins.

To see how galaxies quench their star formation, specifically where in the
galaxies the suppression of star formation first happens and how it progresses
over cosmic time, next we analyze the radial profiles of sSFR. The radial
profiles of sSFR are shown in Figure 8. As we can see from this figure, the
sSFR radial profiles of the majority of our sample galaxies at $z\gtrsim 2.5$
are broadly flat, while they show more diversity in shape at lower redshifts.
At $0.8\lesssim z\lesssim 2.5$, star-forming galaxies in our sample tend to
have a flat or centrally-peaked sSFR$(r)$, while quiescent galaxies tend to
have centrally-suppressed sSFR$(r)$. On the other hand, green-valley galaxies
in our sample seem to have broadly flat radial profiles up to $z\sim 1.0$,
except in the most massive group, where some of them show an sSFR suppression
in their central regions. At lower redshifts, the majority of our sample
galaxies have centrally-suppressed sSFR$(r)$. It is also interesting to see
that the majority of star-forming galaxies at $0.8\lesssim z\lesssim 2.5$ (in
which the cosmic noon epoch is covered), have a centrally-peaked sSFR$(r)$.
This central elevation of sSFR is not observed at higher redshifts.

Figure 8: Similar to Figure 7 but for radial profiles of sSFR. At $z\gtrsim
2.5$, our sample galaxies tend to have a broadly flat sSFR$(r)$, while at the
lower redshifts, they show more diversity in shape. At $0.8\lesssim z\lesssim
2.5$, the majority of star-forming galaxies have an elevation of sSFR in their
central regions, while they tend to have centrally-suppressed sSFR radial
profiles at lower redshifts. On the other hand, quiescent galaxies in our
sample tend to have centrally-suppressed sSFR radial profiles since $z\sim
1.5$.

Next, we analyze the radial profiles of the stellar population age to see how
this quantity varies radially within our sample galaxies and investigate the
underlying stellar population properties causing the diversity in the sSFR
radial profiles. From our spatially resolved SED fitting, we obtain maps of
the mass-weighted ages, which is the average age of stars in a stellar
population as weighted by the stellar mass formed over the course of the star
formation history. The age radial profiles are shown in Figure 9. As can be
seen from this figure, there is a trend of increasing overall age of the
stellar populations in galaxies over cosmic time, as indicated by the
increasing normalization of the radial profiles with decreasing redshift. The
star-forming galaxies that have a centrally-peaked sSFR$(r)$ at $0.8\lesssim
z\lesssim 2.5$ (possibly around the cosmic noon epoch) as shown in Figure 8
are likely in a phase of rapid star formation in their centers (i.e., a
nuclear starburst; e.g., Dekel & Burkert 2014; Zolotov et al. 2015; Tacchella
et al. 2016b; Tadaki et al. 2017), as indicated by the young stellar
populations (age $\lesssim 100$ Myr) in their central regions. At this epoch,
green-valley and quiescent galaxies tend to have radially decreasing age
profiles (i.e., negative gradient). At $0.3<z<0.8$, low-mass galaxies
($\log(M_{*}/M_{\odot})<9.5$) in all stages of star formation have radially
decreasing age radial profiles (i.e., negative gradient). A similar trend
still holds for star-forming and green-valley galaxies in the higher mass
group ($9.5<\log(M_{*}/M_{\odot})<10.5$). On the other hand, quiescent
galaxies at this epoch tend to have overall flat and old stellar populations
across their entire radius, with higher normalization (i.e., older) than that
of star-forming and green-valley galaxies.

Figure 9: Similar to Figure 7 but for radial profiles of the stellar
population age. There is a trend of increasing the overall age of the stellar
populations in galaxies over cosmic time. At $0.8\lesssim z\lesssim 2.5$
(possibly around the cosmic noon epoch), a large fraction of our star-forming
galaxies have a centrally-suppressed age radial profile. Those galaxies have
centrally-peaked sSFR radial profiles as shown in Figure 8. These trends
indicate that they are likely to have an ongoing nuclear starburst event.

###  4.3 Compactness of the Spatial Distributions of Stellar Mass and SFR

The centrally-peaked sSFR$(r)$ of star-forming galaxies at around the cosmic
noon epoch indicates that they are likely undergoing a nuclear starburst that
builds the bulge component, as has also been observed by previous studies
(e.g. Tadaki et al., 2017; Kalita et al., 2022). The centrally-suppressed
sSFR$(r)$ profiles which start to emerge in quiescent galaxies at around the
same epoch can be caused by the cessation of star formation in the center
and/or a matured bulge that has been formed in these galaxies. This trend
provides a hint on how galaxies quench their star formation, which seems to
progress in an inside-to-outside manner (i.e., quenching starts from the
center and then propagates outward). At the same time, this trend may indicate
that galaxies build their central regions first, forming a mature bulge, and
then subsequently assemble their disk through star formation (i.e., inside-out
growth). To further investigate this, next we compare the compactness of the
spatial distributions of $M_{*}$ and SFR by means of the half-mass and half-
SFR radii.

We compare the half-mass radius and the half-SFR radius in Figure 10. The
half-SFR radius is a radius (measured along the elliptical semi-major axis)
that covers half of the total SFR. To compare the distributions of our star-
forming and quiescent galaxies on this diagram, we plot the density contours.
As can be seen from this figure, star-forming galaxies broadly follow the one-
to-one line, whereas quiescent galaxies are in excess above the line. This
means that in quiescent galaxies, the spatial distribution of SFR is more
extended than that of stellar mass, indicating that star formation is ongoing
in the disk and less active in the central region. It is also possible that a
massive bulge might have been formed in the centers, making a more compact
stellar mass distribution. On the other hand, star-forming galaxies are
equally distributed. Some star-forming galaxies have spatially more compact
star formation distribution than the stellar mass (i.e., below the one-to-one
line), which indicates that active star formation happens at their centers. On
the other hand, in the star-forming galaxies that have extended star formation
(i.e., above the one-to-one line), the bulge might have been built and active
star formation is now progressing outward and building the disk.

Figure 10: Comparison between the half-mass radius and half-SFR radius. The
overall symbols in this figure are the same as those in Figure 6. The red and
blue contours show the number densities of the star-forming and quiescent
galaxies. The majority of quiescent galaxies in our sample have spatially more
extended star formation distribution than the stellar mass. Star-forming
galaxies are equally distributed on the diagram. In star-forming galaxies with
compact SFR distribution (i.e., below the one-to-one line), an active star
formation happens in the central region. On the other hand, in star-forming
galaxies with less compact SFR distribution, a bulge might have been built and
active star formation is now progressing outward, building the disk.

It has been known that galaxy size correlates with global $M_{*}$ for galaxies
out to at least $z\sim 3$ (i.e., the size–mass relation; e.g., Shen et al.,
2003; van der Wel et al., 2014; Morishita et al., 2014; Yang et al., 2021).
However, most of the previous studies rely on galaxy half-light radii as a
measure of galaxy size. Since mass-to-light ratios are not constant across a
galaxy’s region but instead have a gradient, the half-light radii are not a
direct probe of the underlying stellar mass profiles. Therefore, it is
expected that the half-mass and half-light radii are different. The difference
in size as probed by light and mass profiles in galaxies out to $z\sim 2.5$
has been investigated by previous studies (e.g., Suess et al., 2019).

Now, we check how half-mass and half-SFR radii correlate with integrated
$M_{*}$ in our sample galaxies. We show the correlations in Figure 11.
Overall, we see that both half-mass and half-SFR radii increase as increasing
integrated $M_{*}$. However, we observe a difference in how star-forming and
quiescent galaxies are distributed in the two correlations. In the top panel,
we can see that star-forming galaxies tend to have a larger half-mass radius
than quiescent galaxies in all masses. On the other hand, as can be seen from
the bottom panel, there is no clear difference between the star-forming and
quiescent galaxies in terms of the half-SFR radius, although star-forming
galaxies tend to have a wider range of half-SFR radius than quiescent
galaxies.

In agreement with the trend we observe here, previous studies using half-light
radii also found that star-forming galaxies are larger than quiescent galaxies
in all $M_{*}$ (van der Wel et al., 2014). However, when the half-mass radius
is used, the normalization of the relation decreases, and the slope as well
becomes shallower (Suess et al., 2019). We do not intend to measure the slope
and normalization of our size–mass relations because of the limited sample.
Ideally, we need a large sample of galaxies and divide them into some redshift
bins. In this work, we combine all galaxies in our sample despite the fact
that they are located in a wide range of redshift.

Figure 11: The half-mass (top panel) and half-SFR (bottom panel) radii as a
function of integrated $M_{*}$. The overall symbols in this figure are the
same as those in Figure 10. The half-mass radius is defined as the radius that
covers half of the total $M_{*}$, while the half-SFR radius is the radius
covering half of the total SFR. Both half-mass and half-SFR radii increase
with increasing $M_{*}$. The $M_{*}$ vs half-mass radius diagram shows that
overall, the star-forming galaxies tend to have a larger size than quiescent
galaxies. However, there is no clear difference between the two groups in
terms of the size of the SFR distribution.

###  4.4 Central Stellar Mass Surface Density

We have observed that quiescent galaxies have a smaller half-mass radius than
star-forming galaxies. The quiescent galaxies at $z\lesssim 2.5$ also have
centrally-suppressed sSFR radial profiles. The compactness of quiescent
galaxies can be caused by the massive bulge that has been formed in their
centers. This massive bulge may also be the reason behind the central
suppression of sSFR. Therefore, here we investigate the stellar mass surface
density in the central $1$ kpc radius ($\Sigma_{*,1\rm{kpc}}$) of our sample
galaxies.

Figure 12 shows a relationship between the integrated $M_{*}$ and
$\Sigma_{*,1\rm{kpc}}$. A tight relationship is evidenced from this figure,
which indicates that $\Sigma_{*,1\rm{kpc}}$ grows hand-in-hand with the
integrated $M_{*}$of the galaxies. We fit the relation involving all galaxies
in our sample with a linear function (in logarithmic scale) using the
Orthogonal Distance Regression (ODR) method and find that the relationship has
a scatter of $0.38$ dex. The best-fit linear function (with a slope of $1.10$
and zero-point of $-2.33$) is shown with the purple dashed line. It is
interesting to see that quiescent galaxies mostly reside in a tight locus at
the top of the relation. This trend indicates that quiescent galaxies tend to
be massive and have a massive central stellar mass density, perhaps associated
with a bulge that has been formed in these galaxies. This massive central
component might in part causes the more compact (i.e., smaller half-mass
radius) size of quiescent galaxies compared to the star-forming galaxies. In
Section 5.1, we will discuss further how this central stellar mass density
evolves with redshift and correlates with star formation in the inner and
outer regions of the galaxies.

A similar relationship has also been observed by previous studies in galaxies
at $z\lesssim 3$ (e.g., Fang et al., 2013; Tacchella et al., 2015; Barro et
al., 2017). Using a larger sample than ours, they also found that quiescent
galaxies occupy a tight locus on top of the overall relationship with all
galaxies. The relation involving only quiescent galaxies has a shallower slope
than that with star-forming galaxies only. The black and red dashed lines are
the best-fit to the $M_{*}$–$\Sigma_{*,1\rm{kpc}}$ relation of $z\sim 0$
passive galaxies reported by Fang et al. (2013) and Tacchella et al. (2015),
respectively, whereas the blue line is the best-fit to the relation of $z\sim
0$ star-forming galaxies from Tacchella et al. (2015). In our result, we also
see a slight bending at the tip of the distribution of quiescent galaxies, as
can be seen from the density contour. It is interesting to see the consistency
between the $M_{*}$–$\Sigma_{*,1\rm{kpc}}$ relation from our study and that
from the literature despite the fact that our sample galaxies cover a wider
redshift range ($0.3<z<6.0$). It suggests that this relation might be
universal and galaxies evolve along this tight relation.

Figure 12: A tight relationship between the integrated $M_{*}$ and stellar
mass surface density within the central $1$ kpc radius
($\Sigma_{*,1\rm{kpc}}$). The overall symbols in this figure are the same as
those in Figure 10. This relation indicates that $\Sigma_{*,1\rm{kpc}}$ grows
hand-in-hand with the integrated $M_{*}$. The quiescent galaxies mostly reside
in the top locus of the relation, indicating that they tend to be massive and
have a massive spheroid in their centers. The purple dashed line is the best-
fit linear function to our sample galaxies, which has a slope of $1.10$, zero-
point of $-2.33$, and scatter of $0.38$ dex. The black and red dashed lines
are the best-fit to the same relation found in $z\sim 0$ passive galaxies
reported by Fang et al. (2013) and Tacchella et al. (2015), respectively,
whereas the blue line is the best-fit to the relation of $z\sim 0$ star-
forming galaxies from Tacchella et al. (2015).

##  5 Discussions

###  5.1 How do Galaxies Grow and Quench Over Cosmic Time?

In this section, we exploit our results to try to infer how galaxies assemble
their structures and eventually cease their star formation activities over the
course of their life. In particular, we are interested to investigate how the
processes of stellar mass buildup and quenching were propagating within the
galaxies. Based on the models of the hierarchical galaxy formation paradigm,
within the framework of the $\Lambda$CDM cosmology, galaxies are predicted to
build their structures over cosmic time in an inside-out manner where the
central component was built first and subsequently the disk structure is
assembled gradually over time (e.g., Cole et al., 2000; van den Bosch, 2002;
Aumer & White, 2013). This assembly process happens through the series of gas
accretions via mergers and filamentary accretion through the cosmic web.

In a gas-rich merger, gas can fall rapidly into the center of the
gravitational potential, causing compaction of gas that eventually triggers a
nuclear starburst event (e.g., Dekel & Burkert, 2014). Other smoother gas
streams can also lead to central gas compaction, e.g., counter-rotating
streams and low angular momentum recycled gas. Based on the zoom-in
cosmological hydrodynamical simulations, this wet compaction event is
predicted to typically occur in $\sim 10^{9.5}M_{\odot}$ galaxies at $z\sim
2-4$ (e.g., Zolotov et al., 2015; Tacchella et al., 2016b, a). During this
event, the star formation at the center is very intense, which converts
massive gas concentration rapidly into stars (i.e., with short depletion
time). This process depletes the gas in galaxies and makes SFR decline.
However, the gas compaction event can occur multiple times in the high
redshift galaxies, causing an up-and-down of their SFRs. This could appear as
an oscillation around the SFMS ridge line. The gas compaction and nuclear
starburst processes can build the massive bulge in the centers of the galaxies
(Tacchella et al., 2016b, a). Subsequently, galaxies can reach full quenching
when the timescale of gas replenishment in the disk is longer than the
timescale of gas depletion by star formation. This is more likely to happen at
lower redshifts when the cosmological gas accretion rate is low, and
especially so in a hot halo above the critical halo mass of $\sim
10^{11.5}M_{\odot}$ (Birnboim & Dekel, 2003; Kereš et al., 2005; Dekel &
Birnboim, 2006). In general, an emerging picture from the above scenarios is
that galaxies build their structures and quench their star formation in an
inside-out manner.

Next, we will exploit our observational results to find indications of the
above model predictions. We note that our current sample is limited, which may
cause some biases in our interpretations. More information from
multiwavelength observations is needed for a more comprehensive study on this,
including maps of the gas mass and their kinematics which can be obtained from
radio and integral field spectroscopy observations. We leave this for future
work.

First, we examine the evolution of the ratio between the half-SFR and half-
mass radii and the ratio between the sSFR in the inner ($\text{sSFR}_{\rm
in}$) and outer ($\text{sSFR}_{\rm out}$) regions of the galaxies. We show the
evolution of these quantities in Figure 13. The $\text{sSFR}_{\rm in}$ and
$\text{sSFR}_{\rm out}$ are defined as the total sSFR inside and outside of
the half-mass radius, respectively. In the study of galaxy evolution, we face
the fact that we observe galaxies at a certain cosmic time. In other words,
the galaxies across the wide redshift that we observe here are not necessarily
connected evolutionarily (i.e., they may not progenitor–descendant pairs),
which makes it difficult to interpret an evolutionary trend from our study.
Previous studies have tried to connect galaxies using a $M_{*}$ growth
function ($M_{*}(z)$) derived from the stellar mass function assuming a
constant number density (e.g., van Dokkum et al., 2010, 2013). To try to
connect galaxies in our sample and infer some evolutionary trends from their
properties, we use $M_{*}(z)$ from van Dokkum et al. (2013) (Eq. 1 therein),
which was designed for tracing the evolution of Milky Way analogs (i.e.,
having $M_{*}=5\times 10^{10}M_{\odot}$ at $z=0$). We bin redshift with a
width of $0.7$. By using the $M_{*}(z)$ function, We then choose galaxies that
fall within $\pm 0.35$ dex in $M_{*}$ at each redshift bin. Our model only
extends up to $z=4.5$ because there are no low-mass galaxies selected beyond
it in our sample. This can be caused by a bias in our sample selection (see
Section 2.2). This toy empirical model is only intended for a reference in
interpreting observational trends. We show the evolutionary trends of this toy
model with black diamond symbols. For comparison, we also show the trends
observed in local spiral galaxies, as inferred from the data analyzed by
Abdurro’uf & Akiyama (2017) and Abdurro’uf et al. (2022a) (personal
communication).

The ratio between the half-SFR and half-mass radii ($R_{e,\rm{SFR}}/R_{e}$)
tells about the relative extent of the SFR distribution compared to the
stellar mass distribution. A ratio $R_{e,\rm{SFR}}/R_{e}$$<1$ implies a more
compact SFR distribution than the stellar mass, whereas
$R_{e,\rm{SFR}}/R_{e}$$>1$ implies a more compact mass distribution than star
formation. $R_{e,\rm{SFR}}/R_{e}$$<<1$ may indicate an ongoing nuclear
starburst, while $R_{e,\rm{SFR}}/R_{e}$$>>1$ indicates that a massive bulge
has been built in the galaxies. From Figure 13, we can see that from an early
epoch up to $z\sim 3.5$, the two ratios are close to unity, indicating a
similar mass doubling time across the galaxy’s region. In our sample, we only
have star-forming and green-valley galaxies at this epoch. Quiescent galaxies
emerge from $z\sim 3$ in our sample and we start to see more dispersion in the
two ratios at this later epoch. At $1.5\lesssim z\lesssim 2.5$ we see a large
fraction of our star-forming galaxies have a compact star formation with
$R_{e,\rm{SFR}}/R_{e}$ as low as $\sim-0.5$ dex and centrally-peaked sSFR with
$\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$ of up to $\sim 3$ dex. In
contrast to this, the majority of green valley and quiescent galaxies at this
epoch have $R_{e,\rm{SFR}}/R_{e}$$>1$ and $\text{sSFR}_{\rm
in}/\text{sSFR}_{\rm out}$$<1$. At a later epoch ($z\lesssim 1.5$), the
majority of quiescent, green valley, and star-forming galaxies have extended
SFR distributions and centrally-suppressed sSFR. The toy model has constant
ratios of $\sim 1$ from the early epoch up to $z\sim 1.5$ after which its sSFR
declines and SFR distribution becomes more extended. There is an indication
that its $\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$ actually increases a
little bit above 1 at $z\sim 1.5$.

The trends observed at $z\sim 0$ indicate that local spiral galaxies have
overall extended SFR distributions and centrally-suppressed sSFR radial
profiles. These trends agree with the scenarios inferred from our results in
the current work and provide a nice extension to our results toward low
redshift, complementing the general picture.

Overall, the above trend agrees with the inside-out growth and quenching
scenarios. In the early cosmic time, galaxies get steady gas accretion for
star formation but they have yet to form a bulge, and the star formation is
likely distributed evenly across their regions. At $z\sim 2$, which coincides
with the peak epoch of the cosmic SFRD and perhaps the cosmic gas accretion
(Madau & Dickinson, 2014), star-forming galaxies in our sample may experience
gas compaction events that later build bulge in their centers. After that,
quenching might have been started in their central regions, but star formation
is still active in the disk that further building the disk. In addition to in-
situ star formation, minor mergers can also contribute to the buildup of
stellar mass in the disk and grow the galaxy size.

Figure 13: Evolution of the ratio between the half-SFR and half-mass radii
(top panel) and the ratio between the sSFR inside and outside of the half-mass
radius (bottom panel). The overall symbols and color-coding are the same as
those in the right panel of Figure 5. The profiles shown with black diamond
symbols represent an expectation from a toy empirical model for the evolution
of the Milky Way analogs. For comparison, we also show the trends observed in
local spiral galaxies from Abdurro’uf & Akiyama (2017) and Abdurro’uf et al.
(2022a). At $z\gtrsim 3.5$, the two ratios are close to unity, indicating a
similar mass doubling time across the galaxy’s region. At $1.5\lesssim
z\lesssim 2.5$, many star-forming galaxies have a compact star formation
($R_{e,\rm{SFR}}/R_{e}$ as low as $\sim-0.5$ dex and $\text{sSFR}_{\rm
in}/\text{sSFR}_{\rm out}$ up to $\sim 3$ dex). The majority of the green
valley and quiescent galaxies at this epoch have $R_{e,\rm{SFR}}/R_{e}$$>1$
and $\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$$<1$. At a later epoch
($z\lesssim 1.5$), the majority of quiescent, green valley, and star-forming
galaxies have extended SFR distributions and centrally-suppressed sSFR.
Overall these trends point toward the inside-out growth and quenching
scenario.

###  5.2 The buildup of the Central Stellar Mass Density Over Cosmic Time

As we have seen in Section 4.4, our sample galaxies exhibit a tight
relationship between the global $M_{*}$ and $\Sigma_{*,1\rm{kpc}}$.
$\Sigma_{*,1\rm{kpc}}$ is a good indicator for quiescent galaxies because they
form a rather distinct sequence at the tip of the overall
$M_{*}$–$\Sigma_{*,1\rm{kpc}}$ relation and have a shallower slope. Here we
discuss the evolution of $\Sigma_{*,1\rm{kpc}}$ with redshift to see how the
central bulge is built over cosmic time in our sample galaxies. As
$\Sigma_{*,1\rm{kpc}}$ develops over time, it is also interesting to analyze
how sSFR at the central 1 kpc evolves following the development of
$\Sigma_{*,1\rm{kpc}}$. The evolution of these two quantities is shown in
Figure 14. As we can see from this figure, $\Sigma_{*,1\rm{kpc}}$ tends to
increase with cosmic time, whereas $\text{sSFR}_{\rm 1kpc}$ declines with
cosmic time. The quiescent galaxies tend to have higher $\Sigma_{*,1\rm{kpc}}$
and lower $\text{sSFR}_{\rm 1kpc}$ in all redshifts. $\text{sSFR}_{\rm 1kpc}$
of quiescent and green-valley galaxies tend to be declined more rapidly than
that of star-forming galaxies. Interestingly, the overall $\text{sSFR}_{\rm
1kpc}$ of our star-forming galaxies does not decline much from $z=6$ up to
$z\sim 1.5$. There is an indication that $\text{sSFR}_{\rm 1kpc}$ of some
star-forming galaxies even increases at $1.5\lesssim z\lesssim 2.5$.

Figure 14: The buildup of the stellar mass density at the central 1 kpc
($\Sigma_{*,1\rm{kpc}}$) over cosmic time (top panel) and the evolution of
sSFR at the central 1 kpc ($\text{sSFR}_{\rm 1kpc}$; bottom panel). The
profiles shown with black diamond symbols represent an expectation from a toy
empirical model for the evolution of the Milky Way analogs, whereas the red
pentagon symbols show a prediction from zoom-in cosmological hydrodynamical
simulations by Tacchella et al. (2016a) for the evolution of massive galaxies
($\log(M_{*}(z=2)/M_{\odot})>10.2$).

The black profiles show the expected evolution of Milky Way analogs based on
our toy model derived in Section 5.1. Its $\Sigma_{*,1\rm{kpc}}$ increases by
$\sim 2$ magnitude over $0.5\lesssim z\lesssim 4.5$, while its
$\text{sSFR}_{\rm 1kpc}$ decreases significantly ($\sim 3.7$ magnitude) over
the same period. This implies that the central SFR within 1 kpc also decreases
with time. At $1.5\lesssim z\lesssim 2.5$, the $\text{sSFR}_{\rm 1kpc}$ of
this model seems to be constant. In addition to our toy model, we also compare
our observational trend with the predictions from the VELA zoom-in
cosmological hydrodynamical simulations (Ceverino et al., 2014; Zolotov et
al., 2015) that was analyzed by Tacchella et al. (2016a). These predictions,
which are shown in red profiles, are obtained by averaging the evolutionary
trends of relatively massive galaxies in the simulations that have
$\log(M_{*}/M_{\odot})=10.2$ at $z=2$. This simulation was run over $1<z<7$.
The $\Sigma_{*,1\rm{kpc}}$ and $\text{sSFR}_{\rm 1kpc}$ predicted from the
cosmological simulation at $z\sim 6$ are in good agreement with our
observations. The evolution of the $\Sigma_{*,1\rm{kpc}}$ and
$\text{sSFR}_{\rm 1kpc}$ from the simulation might be consistent with the
evolution of the progenitors of local massive quiescent galaxies, as implied
from our observations. However, there is an excess of $\Sigma_{*,1\rm{kpc}}$
at $z\sim 0.5$ quiescent galaxies compared to the simulation. Those galaxies
are likely members of the WHL0137$-$08 or MACS0647+70 clusters. In clusters,
galaxies could accrete more mass through mergers which can result in denser
central mass density.

The trends at $z\sim 0$ from Abdurro’uf & Akiyama (2017) and Abdurro’uf et al.
(2022a) provide a good extension for our results in the current work. They
overall agree with the picture of increasing $\Sigma_{*,1\rm{kpc}}$ and
decreasing $\text{sSFR}_{\rm 1kpc}$ with cosmic time. However, we also see a
lower $\Sigma_{*,1\rm{kpc}}$ in the local spiral galaxies than in the
quiescent galaxies at $z\sim 0.5$ that are possibly the cluster members. If we
ignore these galaxies and assume that $\Sigma_{*,1\rm{kpc}}$ trend from the
VELA simulation will evolve to have a similar value as those of the observed
$\Sigma_{*,1\rm{kpc}}$ of local galaxies (which is likely given its shallow
slope at $z\sim 1$), we see a possible saturation of central mass density in
galaxies.

Next, we compare the effects of $\Sigma_{*,1\rm{kpc}}$ and the global $M_{*}$
on $\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$. We show the
$\Sigma_{*,1\rm{kpc}}$–$\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$ and
$M_{*}$–$\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$ relations in Figure 15.
We can see from this figure that $\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$
is correlated with both $\Sigma_{*,1\rm{kpc}}$ and the global $M_{*}$ such
that increasing $\Sigma_{*,1\rm{kpc}}$ and $M_{*}$ corresponds to a steeper
sSFR decline in the central regions. However, the $M_{*}$–$\text{sSFR}_{\rm
in}/\text{sSFR}_{\rm out}$ relation seems to be broader and less significant
compared to $M_{*}$–$\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$, indicating
that $\Sigma_{*,1\rm{kpc}}$ is more influential in driving $\text{sSFR}_{\rm
in}/\text{sSFR}_{\rm out}$ than global $M_{*}$. The majority of those that
have central sSFR suppression are the quiescent and green-valley galaxies,
whereas a significant fraction of star-forming galaxies have broadly flat or
centrally-peaked sSFR radial profile (i.e., negative gradient). This further
suggests that $\Sigma_{*,1\rm{kpc}}$ is a good predictor for quiescent
galaxies, which agree with previous studies (e.g., Fang et al., 2013; Cheung
et al., 2012; Jung et al., 2017; Barro et al., 2017; Whitaker et al., 2017;
Bluck et al., 2022).

Figure 15: The $\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$ ratio as a
function of the central mass density (left panel) and global $M_{*}$(right
panel). The overall symbols in this figure are the same as those in Figure 10.
Galaxies that are massive and have high $\Sigma_{*,1\rm{kpc}}$ tend to have
centrally-suppressed sSFR profiles. This trend is predominantly observed in
quiescent galaxies. The $M_{*}$–$\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$
relation seems to be broader and less significant compared to
$M_{*}$–$\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$, indicating that
$\Sigma_{*,1\rm{kpc}}$ is more influential in driving $\text{sSFR}_{\rm
in}/\text{sSFR}_{\rm out}$ than global $M_{*}$.

Some physical mechanisms for the quenching in galaxies have been proposed,
including the stellar feedback (e.g., Dekel & Silk, 1986; Murray et al.,
2005), AGN feedback (e.g., Di Matteo et al., 2005; Croton et al., 2006; Ciotti
& Ostriker, 2007; Cattaneo et al., 2009), kinematic stabilization of gas in
the disk by massive bulge (e.g., Martig et al., 2009; Genzel et al., 2014),
and a long-term suppression of gas supply due to the virial shock heating by
the dark matter halo when it has reached a threshold mass of $\sim
10^{11.5}M_{\odot}$ (e.g., Birnboim & Dekel, 2003; Dekel & Birnboim, 2006;
Kereš et al., 2009). The inside-out quenching observed in this work seem to
agree with the AGN feedback scenario, as has also been suggested by previous
spatially resolved studies of galaxies at $z\lesssim 2$ (e.g., Bluck et al.,
2022; Nelson et al., 2021). The feedback from AGN can expel or heat up gas,
which then cause the suppression of star formation in at least the central
region of the galaxies. The feedback strength must scale with the mass of the
supermassive black hole (SMBH), which has been known to correlate tightly with
the bulge mass (e.g., Häring & Rix, 2004; Schutte et al., 2019). This is
overall in line with our results which show that quiescent galaxies tend to
have centrally suppressed sSFR and high central mass density. However, more
in-depth study is needed to further investigate this. Ideally, we need multi-
wavelength data sets that allow us to identify AGN and measure the strength of
its feedback, which we have not had in this study.

##  6 Summary and Conclusions

We perform spatially resolved SED fitting on 444 galaxies at $0.3<z<6.0$ in
two clusters (WHL0137$-$08 and MACS0647+70) and a blank field using imaging
data from JWST and HST in up to 13 bands. We use piXedfit throughout the
analysis. This software can simultaneously perform image processing, pixel
binning, and spatially resolved SED fitting. By using the maps of spatially
resolved stellar population properties (on kpc scales) obtained from this
analysis, we investigate how galaxies grow their structures and quench their
star formation activities across cosmic time. Overall, our key results are
summarized in the following:

  1. 1.

The normalization of the stellar mass surface density radial profiles
($\Sigma_{*,1\rm{kpc}}$$(r)$) increases with increasing cosmic time and global
$M_{*}$. At each redshift, quiescent galaxies tend to have higher
$\Sigma_{*,1\rm{kpc}}$ across the entire radius than green-valley and star-
forming galaxies. The sSFR radial profiles (sSFR$(r)$) show more variations
across redshift and global $M_{*}$. The sSFR$(r)$ are broadly flat at
$2.5\lesssim z\lesssim 6.0$ in all galaxies, indicating a similar mass-
doubling time across the entire radius. At $0.8\lesssim z\lesssim 2.5$, less
massive ($\log(M_{*}/M_{\odot})<11.0$) star-forming galaxies have flat or
centrally-peaked sSFR$(r)$, whereas the majority of quiescent galaxies have
centrally-suppressed sSFR$(r)$. At lower redshift ($z<0.8$), almost all
galaxies (regardless of $M_{*}$ and star formation stage) have centrally-
suppressed sSFR$(r)$. The radial profiles of stellar ages show that those
galaxies with centrally-peaked sSFR$(r)$ at $0.8\lesssim z\lesssim 2.5$ have
very young stellar populations in their central regions, indicating an ongoing
nuclear starburst.

  2. 2.

The majority of quiescent galaxies have a larger half-SFR radius than the
half-mass radius, indicating that they have extended spatial distribution of
SFR and compact distribution of stellar mass. In contrast, some star-forming
galaxies, especially at high redshifts, have a half-SFR radius being roughly
similar or smaller to the half-mass radius, whereas those at low redshifts
have a half-SFR radius being larger than the half-mass radius. The half-mass
radius of the star-forming galaxies is on average larger than the quiescent
galaxies in all global $M_{*}$.

  3. 3.

We observe a tight correlation between the global $M_{*}$ and stellar mass
density at the central 1 kpc radius ($\Sigma_{*,1\rm{kpc}}$) with $0.38$ dex,
indicating that galaxies grow their central mass density hand-in-hand with
their global $M_{*}$. The quiescent galaxies reside in a sequence at the tip
of the overall relationship and have a shallower slope. This trend indicates
that $\Sigma_{*,1\rm{kpc}}$ is a good predictor of quenching, where passive
galaxies tend to have higher $\Sigma_{*,1\rm{kpc}}$ and global $M_{*}$. The
shallower slope of $M_{*}$–$\Sigma_{*,1\rm{kpc}}$ in quiescent galaxies
suggests that their central mass density has reached a saturation point.

  4. 4.

We investigate the evolution of the $R_{e,\rm{SFR}}/R_{e}$ and
$\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$ ratios with redshift to try to
understand how galaxies grow their structures and quench their star formations
over cosmic time. We find that the ratios are close to unity from the early
epoch up to $z\sim 3.5$ and the ratios start to deviate from unity since then.
At $1.5\lesssim z\lesssim 2.5$, a fraction of our star-forming sample has a
low $R_{e,\rm{SFR}}/R_{e}$ and high $\text{sSFR}_{\rm in}/\text{sSFR}_{\rm
out}$, indicating that they may be experiencing a nuclear starburst. At the
later epoch, most of our sample galaxies, especially quiescent and green-
valley have a high $R_{e,\rm{SFR}}/R_{e}$ and low $\text{sSFR}_{\rm
in}/\text{sSFR}_{\rm out}$, suggesting that massive bulges might have been
formed in these galaxies and the star formation has been quenched in their
central regions.

  5. 5.

We also investigate the evolution of $\Sigma_{*,1\rm{kpc}}$ and sSFR at the
central 1 kpc ($\text{sSFR}_{\rm 1kpc}$). In general, we see an increasing
$\Sigma_{*,1\rm{kpc}}$ and decreasing $\text{sSFR}_{\rm 1kpc}$ with cosmic
time, indicating the buildup of the central bulge component and the quenching
process in the central region of the galaxies. We also find that quiescent
galaxies tend to have higher $\Sigma_{*,1\rm{kpc}}$ and lower
$\text{sSFR}_{\rm 1kpc}$ than star-forming galaxies in all redshifts.

  6. 6.

Finally, we observe $\Sigma_{*,1\rm{kpc}}$–$\text{sSFR}_{\rm
in}/\text{sSFR}_{\rm out}$ and $M_{*}$–$\text{sSFR}_{\rm in}/\text{sSFR}_{\rm
out}$ relations with a negative slope, indicating that galaxies that are more
massive and have massive $\Sigma_{*,1\rm{kpc}}$ tend to have steeper sSFR
suppression in their centers. The $\Sigma_{*,1\rm{kpc}}$–$\text{sSFR}_{\rm
in}/\text{sSFR}_{\rm out}$ relation seem to be tighter than
$M_{*}$–$\text{sSFR}_{\rm in}/\text{sSFR}_{\rm out}$ indicating that
$\Sigma_{*,1\rm{kpc}}$ is more influential in driving $\text{sSFR}_{\rm
in}/\text{sSFR}_{\rm out}$ than global $M_{*}$. The quiescent galaxies tend to
have higher $\Sigma_{*,1\rm{kpc}}$ and $\text{sSFR}_{\rm in}/\text{sSFR}_{\rm
out}$$<1$, suggesting that the formation of bulge might happen simultaneously
with the quenching of star formation in the central regions.

Our work in this paper demonstrates the great potential of spatially resolved
SED analysis using JWST imaging data. It is interesting to extend this study
with larger sample galaxies taken from various surveys to better understand
the buildup of stellar mass in galaxies and the growth of their structures
over cosmic time. More comprehensive comparisons with zoom-in cosmological
simulations would help to better understand the underlying physics. We will
pursue this in our future work.

## Acknowledgements

We thank the anonymous referee for providing valuable comments that help to
improve this paper. We thank Takahiro Morishita for the useful discussion and
comments. This work is based on observations made with the NASA/ESA/CSA James
Webb Space Telescope (JWST). Some of the data presented in this paper were
obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space
Telescope Science Institute. The specific observations analyzed can be
accessed via http://dx.doi.org/10.17909/d2er-wq71 (catalog
http://dx.doi.org/10.17909/d2er-wq71) and http://dx.doi.org/10.17909/cqfq-5n80
(catalog http://dx.doi.org/10.17909/cqfq-5n80). STScI is operated by the
Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS5–26555. Support to MAST for these data is provided by the NASA
Office of Space Science via grant NAG5–7584 and by other grants and contracts.

A and TH are funded by a grant for JWST-GO-01433 provided by STScI under NASA
contract NAS5-03127. The CosmicDawn Center is funded by the Danish National
Research Foundation (DNRF) under grant #140. PD acknowledges support from the
NWO grant 016.VIDI.189.162 (“ODIN”) and from the European Commission’s and
University of Groningen’s CO-FUND Rosalind Franklin program. RAW acknowledges
support from NASA JWST Interdisciplinary Scientist grants NAG5-12460,
NNX14AN10G and 80NSSC18K0200 from GSFC. AZ and AKM acknowledge support by
grant 2020750 from the United States-Israel Binational Science Foundation
(BSF) and grant 2109066 from the United States National Science Foundation
(NSF), and by the Ministry of Science & Technology, Israel. MO acknowledges
support from JSPS KAKENHI Grant Numbers JP22H01260, JP20H05856, JP20H00181,
JP22K21349. AA acknowledges support from the Swedish Research Council
(Vetenskapsrådet project grants 2021-05559). EV acknowledges financial support
through grants PRIN-MIUR 2017WSCC32, 2020SKSTHZ and the INAF GO Grant 2022
(P.I. E. Vanzella).

##  Appendix A Robustness of the SED Fitting Method: Fitting Tests with Mock
SEDs

To test the robustness of our SED fitting method on this new set of
photometric data, we perform SED-fitting tests using semiempirical mock SEDs,
following a similar procedure as performed in Abdurro’uf et al. (2022a,
Appendix A therein). We draw the parameter values for our mock SEDs from the
measured parameters of real galaxies. In this case, we use the measured
parameters obtained from our fitting to the SEDs within the central effective
radius that were used for determining our photometric redshifts (see Section
3.5). Here, we use 290 galaxies selected randomly from 354 galaxies in the
WHL0137$-$08 (before further exclusion). We prefer to use the parameters of
real galaxies for generating mock SEDs, instead of drawing them randomly
because we can not be sure that the combinations of those random parameters
are physically realistic. Using the set of parameters of 290 galaxies, we
generate mock SEDs using the same modeling setup as in our main analysis. We
generate two sets of mock SEDs. The first one with 12 filters of JWST and HST,
the same set of filters as available for the WHL0137$-$08 cluster. The second
set with 8 JWST filters, excluding HST filters. We then inject Gaussian noises
assuming S/N of 20 in all filters. After that, we fit the mock SEDs using the
same method as we used in the main analysis of this paper. For simplicity,
here we fix redshift.

We present the results in Figure 16, which shows the comparisons between the
best-fit parameters derived from SED fitting and the true values from the mock
SEDs. Histograms in the insets show ratios between the best-fit parameters and
the true values. Results of the SED fitting with two sets of photometry are
shown with different symbols and the histograms are shown with different
colors. The data points are color-coded based on their redshifts.

Overall, our SED fitting can recover the true parameters reasonably well.
Stellar mass and mass-weighted age are recovered very well in both 8-band and
12-band SED fitting with small offset ($\lesssim 0.06$ dex) and small standard
deviation ($\lesssim 0.3$ dex). It is interesting to see that mass-weighed age
is well recovered here. This may be due to the fact that our photometry covers
the Balmer break and we have good photometry in rest-frame NIR from JWST,
which also provides a good constraint for $M_{*}$ (see Appendix B). The
stellar metallicity ($Z$) and dust attenuation in the diffuse ISM ($A_{V,2}$)
are also recovered well with an offset of $\lesssim 0.07$ dex and a standard
deviation of $\lesssim 0.5$ dex although they look to be more scattered due to
its small dynamical range. The SFR is more difficult to be recovered for
passive galaxies ($\log(\text{SFR})\lesssim-1$) than for star-forming ones.
For the whole sample, the SED fitting with 12 bands gives better SFR estimates
(a small offset of $0.07$ dex and scatter of $0.64$ dex) than with only JWST
bands (offset of $0.16$ dex and scatter of $0.77$ dex).

Figure 16: Comparisons between the best-fit parameters derived from SED-
fitting tests using mock SEDs and the ground truth. In this test, two sets of
photometry are used, one with combined 12 filters of HST and JWST and the
other with 8 filters of JWSTṪhe results of those two SED fitting tests are
shown with different symbols (circles and triangles). Histograms in the insets
show ratios between the best-fit parameters and the true values from the mock
SEDs. The color-coding of the data points represents redshift.

##  Appendix B The Age–Dust–Metallicity Degeneracy

The addition of JWST NIRCam data extends the wavelength coverage up to roughly
the rest-frame NIR for our sample galaxies. This sufficiently wide wavelength
coverage has the potential to break the well-known age–dust–metallicity
degeneracy in SED fitting, which is very important for the analysis in this
paper. In particular, it is crucial to be able to determine if the reddening
observed in the central region of some galaxies in our sample is due to aging
(i.e., quiescence) or dust attenuation. To check if our data set provides a
sufficient constraint for resolving this degeneracy, we examine model color-
color diagrams among the NIRCam filters. We base our analysis on the observed
frame, instead of the rest-frame, and explore different sets of filters for
different redshifts. To generate model SEDs, we use an overall similar setting
as that used in the main analysis of this paper and assume a double power-law
SFH with $\alpha=3.0$, $\beta=0.5$, and $\tau=1.5$ Gyr, and
$M_{*}=10^{10.5}M_{\odot}$. We then generate model SEDs in grids of age,
$A_{V,2}$, and metallicity. For simplicity, we assume $A_{V,1}=1.5\times
A_{V,2}$.

Figure 17 shows the color-color diagrams of models at $z=2.5$ and $3.5$.
Different colors represent different ages, whereas increasing marker size
represents increasing metallicity (for the square symbol in the first column)
and $A_{V,2}$ (for the circle symbol in the second column). For the two plots
in the first column, we fix $A_{V,2}=0.0$, while for the two plots in the
second column, we fix $Z=Z_{\odot}$. F090W$-$F200W vs. F200W$-$F444W diagram
at the two redshifts seems to be able to distinguish the reddening effect by
age and metallicity in such a way that both effects are almost orthogonal with
each other. An increasing $Z$ at a fixed age corresponds to reddening in
F200W$-$F444W and roughly constant F090W$-$F200W (i.e., vertical shift on the
diagram). On the other hand, an increasing age tends to make reddening in the
two colors (i.e., a diagonal shift on the diagram) for galaxies with low $Z$,
while it corresponds to a roughly horizontal shift for galaxies with high $Z$.

In the second column, we show a relative effect of the dust attenuation and
aging on F090W$-$F277W vs. F277W$-$F444W (for $z=2.5$) and F150W$-$F356W vs.
F356W$-$F444W (for $z=3.5$) diagrams. The two effects are distinguishable in
these two diagrams with dust attenuation seems to make a diagonal shift,
whereas an aging effect is roughly orthogonal to it. For an illustration, in
the rightmost panel, we show examples of the model SEDs of an old less-dusty
(brown color; age$1.5$ Gyr, $\log(Z/Z_{\odot})=-0.5$, $A_{V,2}=0.1$ mag),
young metal-rich (purple; age$0.1$ Gyr, $\log(Z/Z_{\odot})=0.2$, $A_{V,2}=0.1$
mag), and young dusty (green; age$0.1$ Gyr, $\log(Z/Z_{\odot})=-0.5$,
$A_{V,2}=2.0$ mag) galaxies $z=3.5$. We normalize the SEDs by dividing them
with the F200W flux. We also show the transmission curves of the filters that
are used in the color-color diagrams for $z=3.5$ in the left two columns. As
we can see from this figure, the three model SEDs are distinguishable with
NIRCam photometry. The reddening due to the dust attenuation is easily
recognizable in the rest-frame UV to NIR, while due to the metallicity is not
easily recognizable in the rest-frame UV colors but it is detectable in the
rest-frame around the Balmer break ($\sim 4000$ Å) and NIR.

Figure 17: The patterns of model stellar populations on the observed-frame
color-color diagrams involving JWST NIRCam filters. We plot the patterns of
models at $z=2.5$ and $z=3.5$ on the color-color diagrams to check how good
the NIRCam photometry is in resolving the degeneracies among age, dust
attenuation, and metallicity in high redshifts. Different colors represent
different ages, whereas an increasing marker size represents an increasing $Z$
(for the square symbol) and increasing $A_{V,2}$ (for the circle symbol). The
rightmost panel shows examples of the model SEDs of old less-dusty (brown
color), young metal-rich (purple), and young dusty (green) galaxies at
$z=3.5$. The reddening effects by aging and dust are almost orthogonal with
each other (i.e., distinguishable) on the F090W$-$F200W vs. F200W$-$F444W
diagram. The reddening effects by aging and dust attenuation are
distinguishable on the F090W$-$F277W vs. F277W$-$F444W for $z=2.5$ and the
F150W$-$F356W vs. F356W$-$F444W for $z=3.5$.

##  Appendix C Construction of empirical Point Spread Functions and
Convolution Kernels

We generate the empirical PSFs of the HST ACS and JWST NIRcam filters in the
WHL0137$-$08 and MACS0647+70 clusters by stacking images of bright isolated
stars in those fields. For this, we use photutils package (Bradley et al.,
2022a). We show the encircled energy of the empirical PSFs in the left column
of Figure 18. After generating the PSFs, we then construct the convolution
kernels that can be used for PSF matching. As described in Section 3.2, we
perform PSF matching to homogenize the PSF sizes of our imaging data to match
the PSF size of F444W, which is the largest among the filters used in our
work. We also use photutils for generating the kernels. To check the
reliability of our kernels and PSF matching, we convolve the PSF images of the
filters other than F444W with the kernels and compare the encircled energy of
the convolved-PSFs to that of F444W filter. We show this comparison in the
right column of Figure 18. As we can see from this figure, the convolved PSFs
have similar encircled energy (with a small deviation $<0.1$ dex around a
radius of $0\farcs 1$), indicating the robustness of our PSF matching process.

Figure 18: Left column: Encircled energy of the empirical PSFs of the HST ACS
and JWST NIRCam filters in WHL0137$-$08 and MACS0647+70 clusters. The PSFs are
generated by stacking the images of bright isolated stars in the fields. Right
column: Encircled energy of the PSFs after being convolved with kernels. The
convolved PSFs have similar encircled energy, indicating the robustness of our
PSF matching process.

## References

  * Abdurro’uf & Akiyama (2017) Abdurro’uf, & Akiyama, M. 2017, MNRAS, 469, 2806, doi: 10.1093/mnras/stx936
  * Abdurro’uf & Akiyama (2018) —. 2018, MNRAS, 479, 5083, doi: 10.1093/mnras/sty1771
  * Abdurro’uf et al. (2022a) Abdurro’uf, Lin, Y.-T., Hirashita, H., et al. 2022a, ApJ, 926, 81, doi: 10.3847/1538-4357/ac439a
  * Abdurro’uf et al. (2022b) —. 2022b, ApJ, 935, 98, doi: 10.3847/1538-4357/ac7da4
  * Abdurro’uf et al. (2021) Abdurro’uf, Lin, Y.-T., Wu, P.-F., & Akiyama, M. 2021, ApJS, 254, 15, doi: 10.3847/1538-4365/abebe2
  * Abdurro’uf et al. (2022c) —. 2022c, piXedfit: Analyze spatially resolved SEDs of galaxies, Astrophysics Source Code Library, record ascl:2207.033.  http://ascl.net/2207.033
  * 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
  * Aumer & White (2013) Aumer, M., & White, S. D. M. 2013, MNRAS, 428, 1055, doi: 10.1093/mnras/sts083
  * Baker et al. (2022) Baker, W. M., Maiolino, R., Bluck, A. F. L., et al. 2022, MNRAS, 510, 3622, doi: 10.1093/mnras/stab3672
  * Barbary (2016) Barbary, K. 2016, Journal of Open Source Software, 1, 58, doi: 10.21105/joss.00058
  * Barro et al. (2017) Barro, G., Faber, S. M., Koo, D. C., et al. 2017, ApJ, 840, 47, doi: 10.3847/1538-4357/aa6b05
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Birnboim & Dekel (2003) Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349, doi: 10.1046/j.1365-8711.2003.06955.x
  * Bluck et al. (2022) Bluck, A. F. L., Maiolino, R., Brownson, S., et al. 2022, A&A, 659, A160, doi: 10.1051/0004-6361/202142643
  * Bluck et al. (2020) Bluck, A. F. L., Maiolino, R., Piotrowska, J. M., et al. 2020, MNRAS, 499, 230, doi: 10.1093/mnras/staa2806
  * Bouché et al. (2010) Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001, doi: 10.1088/0004-637X/718/2/1001
  * Boyer et al. (2022) Boyer, M. L., Anderson, J., Gennaro, M., et al. 2022, arXiv e-prints, arXiv:2209.03348.  https://arxiv.org/abs/2209.03348
  * Bradley et al. (2022a) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022a, astropy/photutils: 1.5.0, 1.5.0, Zenodo, Zenodo, doi: 10.5281/zenodo.6825092
  * Bradley et al. (2022b) Bradley, L. D., Coe, D., Brammer, G., et al. 2022b, arXiv e-prints, arXiv:2210.01777.  https://arxiv.org/abs/2210.01777
  * Brammer (2022) Brammer, G. 2022, Preliminary updates to the NIRCam photometric calibration, Zenodo, doi: 10.5281/zenodo.7143382
  * Brammer et al. (2022) Brammer, G., Strait, V., Matharu, J., & Momcheva, I. 2022, grizli, 1.5.0, Zenodo, doi: 10.5281/zenodo.6672538
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151, doi: 10.1111/j.1365-2966.2004.07881.x
  * Broadhurst et al. (2005) Broadhurst, T., Benítez, N., Coe, D., et al. 2005, ApJ, 621, 53, doi: 10.1086/426494
  * Bundy et al. (2015) Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ, 798, 7, doi: 10.1088/0004-637X/798/1/7
  * Cano-Díaz et al. (2016) Cano-Díaz, M., Sánchez, S. F., Zibetti, S., et al. 2016, ApJ, 821, L26, doi: 10.3847/2041-8205/821/2/L26
  * Casey et al. (2018) Casey, C. M., Zavala, J. A., Spilker, J., et al. 2018, ApJ, 862, 77, doi: 10.3847/1538-4357/aac82d
  * Cattaneo et al. (2009) Cattaneo, A., Faber, S. M., Binney, J., et al. 2009, Nature, 460, 213, doi: 10.1038/nature08135
  * Ceverino et al. (2014) Ceverino, D., Klypin, A., Klimek, E. S., et al. 2014, MNRAS, 442, 1545, doi: 10.1093/mnras/stu956
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Chan et al. (2017) Chan, B. M. Y., Broadhurst, T., Lim, J., et al. 2017, ApJ, 835, 44, doi: 10.3847/1538-4357/835/1/44
  * Charlot & Fall (2000) Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718, doi: 10.1086/309250
  * Chen et al. (2022) Chen, C.-C., Gao, Z.-K., Hsu, Q.-N., et al. 2022, ApJ, 939, L7, doi: 10.3847/2041-8213/ac98c6
  * Cheung et al. (2012) Cheung, E., Faber, S. M., Koo, D. C., et al. 2012, ApJ, 760, 131, doi: 10.1088/0004-637X/760/2/131
  * Ciotti & Ostriker (2007) Ciotti, L., & Ostriker, J. P. 2007, ApJ, 665, 1038, doi: 10.1086/519833
  * Coe et al. (2008) Coe, D., Fuselier, E., Benítez, N., et al. 2008, ApJ, 681, 814, doi: 10.1086/588250
  * Coe et al. (2012) Coe, D., Umetsu, K., Zitrin, A., et al. 2012, ApJ, 757, 22, doi: 10.1088/0004-637X/757/1/22
  * Coe et al. (2013) Coe, D., Zitrin, A., Carrasco, M., et al. 2013, ApJ, 762, 32, doi: 10.1088/0004-637X/762/1/32
  * Coe et al. (2019) Coe, D., Salmon, B., Bradač, M., et al. 2019, ApJ, 884, 85, doi: 10.3847/1538-4357/ab412b
  * Cole et al. (2000) Cole, S., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2000, MNRAS, 319, 168, doi: 10.1046/j.1365-8711.2000.03879.x
  * Conroy & Gunn (2010) Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833, doi: 10.1088/0004-637X/712/2/833
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486
  * Croton et al. (2006) Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11, doi: 10.1111/j.1365-2966.2005.09675.x
  * Daddi et al. (2007) Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156, doi: 10.1086/521818
  * Daddi et al. (2010) Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686, doi: 10.1088/0004-637X/713/1/686
  * Dayal & Ferrara (2018) Dayal, P., & Ferrara, A. 2018, Phys. Rep., 780, 1, doi: 10.1016/j.physrep.2018.10.002
  * Dekel & Birnboim (2006) Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2, doi: 10.1111/j.1365-2966.2006.10145.x
  * Dekel & Burkert (2014) Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870, doi: 10.1093/mnras/stt2331
  * Dekel & Silk (1986) Dekel, A., & Silk, J. 1986, ApJ, 303, 39, doi: 10.1086/164050
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604, doi: 10.1038/nature03335
  * Diego et al. (2005) Diego, J. M., Protopapas, P., Sandvik, H. B., & Tegmark, M. 2005, MNRAS, 360, 477, doi: 10.1111/j.1365-2966.2005.09021.x
  * Diego et al. (2007) Diego, J. M., Tegmark, M., Protopapas, P., & Sandvik, H. B. 2007, MNRAS, 375, 958, doi: 10.1111/j.1365-2966.2007.11380.x
  * Dimauro et al. (2022) Dimauro, P., Daddi, E., Shankar, F., et al. 2022, MNRAS, 513, 256, doi: 10.1093/mnras/stac884
  * Ebeling et al. (2007) Ebeling, H., Barrett, E., Donovan, D., et al. 2007, ApJ, 661, L33, doi: 10.1086/518603
  * Elbaz et al. (2007) Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33, doi: 10.1051/0004-6361:20077525
  * Ellison et al. (2018) Ellison, S. L., Sánchez, S. F., Ibarra-Medel, H., et al. 2018, MNRAS, 474, 2039, doi: 10.1093/mnras/stx2882
  * Enia et al. (2020) Enia, A., Rodighiero, G., Morselli, L., et al. 2020, MNRAS, 493, 4107, doi: 10.1093/mnras/staa433
  * Falcón-Barroso et al. (2011) Falcón-Barroso, J., Sánchez-Blázquez, P., Vazdekis, A., et al. 2011, A&A, 532, A95, doi: 10.1051/0004-6361/201116842
  * Fang et al. (2013) Fang, J. J., Faber, S. M., Koo, D. C., & Dekel, A. 2013, ApJ, 776, 63, doi: 10.1088/0004-637X/776/1/63
  * Ferland et al. (1998) Ferland, G. J., Korista, K. T., Verner, D. A., et al. 1998, PASP, 110, 761, doi: 10.1086/316190
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mexicana Astron. Astrofis., 49, 137.  https://arxiv.org/abs/1302.4485
  * Ferreira et al. (2022) Ferreira, L., Adams, N., Conselice, C. J., et al. 2022, ApJ, 938, L2, doi: 10.3847/2041-8213/ac947c
  * Finkelstein et al. (2022) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2022, arXiv e-prints, arXiv:2211.05792.  https://arxiv.org/abs/2211.05792
  * Förster Schreiber et al. (2018) Förster Schreiber, N. M., Renzini, A., Mancini, C., et al. 2018, ApJS, 238, 21, doi: 10.3847/1538-4365/aadd49
  * Fudamoto et al. (2021) Fudamoto, Y., Oesch, P. A., Schouws, S., et al. 2021, Nature, 597, 489, doi: 10.1038/s41586-021-03846-z
  * 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
  * Genzel et al. (2010) Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091, doi: 10.1111/j.1365-2966.2010.16969.x
  * Genzel et al. (2014) Genzel, R., Förster Schreiber, N. M., Lang, P., et al. 2014, ApJ, 785, 75, doi: 10.1088/0004-637X/785/1/75
  * Giménez-Arteaga et al. (2022) Giménez-Arteaga, C., Oesch, P. A., Brammer, G. B., et al. 2022, arXiv e-prints, arXiv:2212.08670.  https://arxiv.org/abs/2212.08670
  * Girardi et al. (2000) Girardi, L., Bressan, A., Bertelli, G., & Chiosi, C. 2000, A&AS, 141, 371, doi: 10.1051/aas:2000126
  * Häring & Rix (2004) Häring, N., & Rix, H.-W. 2004, ApJ, 604, L89, doi: 10.1086/383567
  * Hoffmann et al. (2021) Hoffmann, S. L., Mack, J., Avila, R., et al. 2021, in American Astronomical Society Meeting Abstracts, Vol. 53, American Astronomical Society Meeting Abstracts, 216.02
  * Hsiao et al. (2022) Hsiao, T. Y.-Y., Coe, D., Abdurro’uf, et al. 2022, arXiv e-prints, arXiv:2210.14123.  https://arxiv.org/abs/2210.14123
  * Hsieh et al. (2017) Hsieh, B. C., Lin, L., Lin, J. H., et al. 2017, ApJ, 851, L24, doi: 10.3847/2041-8213/aa9d80
  * Inoue et al. (2014) Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805, doi: 10.1093/mnras/stu936
  * Iyer et al. (2018) Iyer, K., Gawiser, E., Davé, R., et al. 2018, ApJ, 866, 120, doi: 10.3847/1538-4357/aae0fa
  * Jullo & Kneib (2009) Jullo, E., & Kneib, J. P. 2009, MNRAS, 395, 1319, doi: 10.1111/j.1365-2966.2009.14654.x
  * Jullo et al. (2007) Jullo, E., Kneib, J. P., Limousin, M., et al. 2007, New Journal of Physics, 9, 447, doi: 10.1088/1367-2630/9/12/447
  * Jung et al. (2017) Jung, I., Finkelstein, S. L., Song, M., et al. 2017, ApJ, 834, 81, doi: 10.3847/1538-4357/834/1/81
  * Kalita et al. (2022) Kalita, B. S., Daddi, E., Bournaud, F., et al. 2022, A&A, 666, A44, doi: 10.1051/0004-6361/202243100
  * Kartaltepe et al. (2022) Kartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2022, arXiv e-prints, arXiv:2210.14713.  https://arxiv.org/abs/2210.14713
  * Kereš et al. (2009) Kereš, D., Katz, N., Fardal, M., Davé, R., & Weinberg, D. H. 2009, MNRAS, 395, 160, doi: 10.1111/j.1365-2966.2009.14541.x
  * Kereš et al. (2005) Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2, doi: 10.1111/j.1365-2966.2005.09451.x
  * Kneib et al. (1993) Kneib, J. P., Mellier, Y., Fort, B., & Mathez, G. 1993, A&A, 273, 367
  * Koekemoer et al. (2003) Koekemoer, A. M., Fruchter, A. S., Hook, R. N., & Hack, W. 2003, in HST Calibration Workshop : Hubble after the Installation of the ACS and the NICMOS Cooling System, 337
  * Leja et al. (2022) Leja, J., Speagle, J. S., Ting, Y.-S., et al. 2022, ApJ, 936, 165, doi: 10.3847/1538-4357/ac887d
  * Leslie et al. (2020) Leslie, S. K., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58, doi: 10.3847/1538-4357/aba044
  * Lin et al. (2019) Lin, L., Pan, H.-A., Ellison, S. L., et al. 2019, ApJ, 884, L33, doi: 10.3847/2041-8213/ab4815
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Marigo & Girardi (2007) Marigo, P., & Girardi, L. 2007, A&A, 469, 239, doi: 10.1051/0004-6361:20066772
  * Marigo et al. (2008) Marigo, P., Girardi, L., Bressan, A., et al. 2008, A&A, 482, 883, doi: 10.1051/0004-6361:20078467
  * Martig et al. (2009) Martig, M., Bournaud, F., Teyssier, R., & Dekel, A. 2009, ApJ, 707, 250, doi: 10.1088/0004-637X/707/1/250
  * McGreer et al. (2015) McGreer, I. D., Mesinger, A., & D’Odorico, V. 2015, MNRAS, 447, 499, doi: 10.1093/mnras/stu2449
  * Meena et al. (2022) Meena, A. K., Zitrin, A., Jiménez-Teja, Y., et al. 2022, arXiv e-prints, arXiv:2211.13334.  https://arxiv.org/abs/2211.13334
  * Morishita et al. (2014) Morishita, T., Ichikawa, T., & Kajisawa, M. 2014, ApJ, 785, 18, doi: 10.1088/0004-637X/785/1/18
  * Morishita et al. (2015) Morishita, T., Ichikawa, T., Noguchi, M., et al. 2015, ApJ, 805, 34, doi: 10.1088/0004-637X/805/1/34
  * Morrissey et al. (2007) Morrissey, P., Conrow, T., Barlow, T. A., et al. 2007, ApJS, 173, 682, doi: 10.1086/520512
  * Morselli et al. (2020) Morselli, L., Rodighiero, G., Enia, A., et al. 2020, MNRAS, 496, 4606, doi: 10.1093/mnras/staa1811
  * Murray et al. (2005) Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569, doi: 10.1086/426067
  * Nardiello et al. (2022) Nardiello, D., Bedin, L. R., Burgasser, A., et al. 2022, arXiv e-prints, arXiv:2209.06547.  https://arxiv.org/abs/2209.06547
  * Nelson et al. (2012) Nelson, E. J., van Dokkum, P. G., Brammer, G., et al. 2012, ApJ, 747, L28, doi: 10.1088/2041-8205/747/2/L28
  * Nelson et al. (2016) Nelson, E. J., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2016, ApJ, 828, 27, doi: 10.3847/0004-637X/828/1/27
  * Nelson et al. (2021) Nelson, E. J., Tacchella, S., Diemer, B., et al. 2021, MNRAS, 508, 219, doi: 10.1093/mnras/stab2131
  * Noeske et al. (2007) Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43, doi: 10.1086/517926
  * Oguri (2010) Oguri, M. 2010, PASJ, 62, 1017, doi: 10.1093/pasj/62.4.1017
  * Pascale et al. (2022) Pascale, M., Frye, B. L., Diego, J., et al. 2022, ApJ, 938, L6, doi: 10.3847/2041-8213/ac9316
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266, doi: 10.1086/340952
  * Pérez-González et al. (2022) Pérez-González, P. G., Barro, G., Annunziatella, M., et al. 2022, arXiv e-prints, arXiv:2211.00045.  https://arxiv.org/abs/2211.00045
  * Popesso et al. (2023) Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526, doi: 10.1093/mnras/stac3214
  * Postman et al. (2012) Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25, doi: 10.1088/0067-0049/199/2/25
  * Rieke et al. (2022) Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2022, arXiv e-prints, arXiv:2212.12069.  https://arxiv.org/abs/2212.12069
  * Rigby et al. (2022) Rigby, J., Perrin, M., McElwain, M., et al. 2022, arXiv e-prints, arXiv:2207.05632.  https://arxiv.org/abs/2207.05632
  * Roberts-Borsani et al. (2022) Roberts-Borsani, G., Treu, T., Chen, W., et al. 2022, arXiv e-prints, arXiv:2210.15639.  https://arxiv.org/abs/2210.15639
  * Salmon et al. (2015) Salmon, B., Papovich, C., Finkelstein, S. L., et al. 2015, ApJ, 799, 183, doi: 10.1088/0004-637X/799/2/183
  * Sánchez (2020) Sánchez, S. F. 2020, ARA&A, 58, 99, doi: 10.1146/annurev-astro-012120-013326
  * Sánchez et al. (2013) Sánchez, S. F., Rosales-Ortega, F. F., Jungwiert, B., et al. 2013, A&A, 554, A58, doi: 10.1051/0004-6361/201220669
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, MNRAS, 371, 703, doi: 10.1111/j.1365-2966.2006.10699.x
  * Santini et al. (2017) Santini, P., Fontana, A., Castellano, M., et al. 2017, ApJ, 847, 76, doi: 10.3847/1538-4357/aa8874
  * Schutte et al. (2019) Schutte, Z., Reines, A. E., & Greene, J. E. 2019, ApJ, 887, 245, doi: 10.3847/1538-4357/ab35dd
  * Shen et al. (2003) Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978, doi: 10.1046/j.1365-8711.2003.06740.x
  * Skelton et al. (2014) Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24, doi: 10.1088/0067-0049/214/2/24
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15, doi: 10.1088/0067-0049/214/2/15
  * Suess et al. (2019) Suess, K. A., Kriek, M., Price, S. H., & Barro, G. 2019, ApJ, 877, 103, doi: 10.3847/1538-4357/ab1bda
  * Tacchella et al. (2016a) Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016a, MNRAS, 458, 242, doi: 10.1093/mnras/stw303
  * Tacchella et al. (2016b) —. 2016b, MNRAS, 457, 2790, doi: 10.1093/mnras/stw131
  * Tacchella et al. (2020) Tacchella, S., Forbes, J. C., & Caplar, N. 2020, MNRAS, 497, 698, doi: 10.1093/mnras/staa1838
  * Tacchella et al. (2015) Tacchella, S., Carollo, C. M., Renzini, A., et al. 2015, Science, 348, 314, doi: 10.1126/science.1261094
  * Tacchella et al. (2018) Tacchella, S., Carollo, C. M., Förster Schreiber, N. M., et al. 2018, ApJ, 859, 56, doi: 10.3847/1538-4357/aabf8b
  * Tadaki et al. (2017) Tadaki, K.-i., Genzel, R., Kodama, T., et al. 2017, ApJ, 834, 135, doi: 10.3847/1538-4357/834/2/135
  * Tomczak et al. (2016) Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2016, ApJ, 817, 118, doi: 10.3847/0004-637X/817/2/118
  * Treu et al. (2013) Treu, T., Schmidt, K. B., Trenti, M., Bradley, L. D., & Stiavelli, M. 2013, ApJ, 775, L29, doi: 10.1088/2041-8205/775/1/L29
  * van den Bosch (2002) van den Bosch, F. C. 2002, MNRAS, 332, 456, doi: 10.1046/j.1365-8711.2002.05328.x
  * van der Wel et al. (2014) van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28, doi: 10.1088/0004-637X/788/1/28
  * van Dokkum et al. (2010) van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018, doi: 10.1088/0004-637X/709/2/1018
  * van Dokkum et al. (2013) van Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJ, 771, L35, doi: 10.1088/2041-8205/771/2/L35
  * van Dokkum et al. (2015) van Dokkum, P. G., Nelson, E. J., Franx, M., et al. 2015, ApJ, 813, 23, doi: 10.1088/0004-637X/813/1/23
  * Vanzella et al. (2022) Vanzella, E., Claeyssens, A., Welch, B., et al. 2022, arXiv e-prints, arXiv:2211.09839.  https://arxiv.org/abs/2211.09839
  * Welch et al. (2022a) Welch, B., Coe, D., Diego, J. M., et al. 2022a, Nature, 603, 815, doi: 10.1038/s41586-022-04449-y
  * Welch et al. (2022b) Welch, B., Coe, D., Zackrisson, E., et al. 2022b, arXiv e-prints, arXiv:2208.09007.  https://arxiv.org/abs/2208.09007
  * Wen & Han (2015) Wen, Z. L., & Han, J. L. 2015, ApJ, 807, 178, doi: 10.1088/0004-637X/807/2/178
  * Wen et al. (2012) Wen, Z. L., Han, J. L., & Liu, F. S. 2012, ApJS, 199, 34, doi: 10.1088/0067-0049/199/2/34
  * Whitaker et al. (2012) Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29, doi: 10.1088/2041-8205/754/2/L29
  * Whitaker et al. (2014) Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104, doi: 10.1088/0004-637X/795/2/104
  * Whitaker et al. (2017) Whitaker, K. E., Bezanson, R., van Dokkum, P. G., et al. 2017, ApJ, 838, 19, doi: 10.3847/1538-4357/aa6258
  * Williams et al. (2022) Williams, H., Kelly, P. L., Chen, W., et al. 2022, arXiv e-prints, arXiv:2210.15699.  https://arxiv.org/abs/2210.15699
  * Wuyts et al. (2013) Wuyts, S., Förster Schreiber, N. M., Nelson, E. J., et al. 2013, ApJ, 779, 135, doi: 10.1088/0004-637X/779/2/135
  * Yang et al. (2021) Yang, L., Roberts-Borsani, G., Treu, T., et al. 2021, MNRAS, 501, 1028, doi: 10.1093/mnras/staa3713
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513
  * Zitrin et al. (2011) Zitrin, A., Broadhurst, T., Barkana, R., Rephaeli, Y., & Benítez, N. 2011, MNRAS, 410, 1939, doi: 10.1111/j.1365-2966.2010.17574.x
  * Zitrin et al. (2009) Zitrin, A., Broadhurst, T., Umetsu, K., et al. 2009, MNRAS, 396, 1985, doi: 10.1111/j.1365-2966.2009.14899.x
  * Zitrin et al. (2015) Zitrin, A., Fabris, A., Merten, J., et al. 2015, ApJ, 801, 44, doi: 10.1088/0004-637X/801/1/44
  * Zolotov et al. (2015) Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327, doi: 10.1093/mnras/stv740