# JADES + JEMS: A Detailed Look at the Buildup of Central Stellar Cores and Suppression of Star Formation in Galaxies at Redshifts $3 < z < 4.5$

Zhiyuan Ji<fnref n="1" />, Christina C. Williams<fnref n="1" />,<fnref n="2" />, Sandro Tacchella<fnref n="3" />,<fnref n="4" />, Katherine A. Suess<fnref n="5" />,<fnref n="6" />, William M. Baker<fnref n="3" />,<fnref n="4" />, Stacey Alberts<fnref n="1" />, Andrew J. Bunker<fnref n="7" />, Benjamin D. Johnson<fnref n="8" />, Brant Robertson<fnref n="5" />, Fengwu Sun<fnref n="1" />, Daniel J. Eisenstein<fnref n="8" />, Marcia Rieke<fnref n="1" />, Michael V. Maseda<fnref n="9" />, Kevin Hainline<fnref n="1" />, Ryan Hausen<fnref n="10" />, George Rieke<fnref n="1" />, Christopher N. A. Willmer<fnref n="1" />, Eiichi Egami<fnref n="1" />, Irene Shivaei<fnref n="1" />, Stefano Carniani<fnref n="11" />, Stephane Charlot<fnref n="12" />, Jacopo Chevallard<fnref n="7" />, Emma Curtis-Lake<fnref n="13" />, Tobias J. Looser<fnref n="3" />,<fnref n="4" />, Roberto Maiolino<fnref n="3" />,<fnref n="4" />,<fnref n="14" />, Chris Willott<fnref n="15" />, Zuyi Chen<fnref n="1" />, Jakob M. Helton<fnref n="1" />, Jianwei Lyu<fnref n="1" />, Erica Nelson<fnref n="16" />, Rachana Bhatawdekar<fnref n="17" />,<fnref n="18" />, Kristan Boyett<fnref n="19" />,<fnref n="20" />, and Lester Sandles<fnref n="3" />,<fnref n="4" />

## Abstract

We present a spatially resolved study of stellar populations in six galaxies with stellar masses $M_* \sim 10^{10} M_\odot$ at $z \sim 3.7$ using 14-filter James Webb Space Telescope (JWST)/NIRCam imaging from the JADES and JEMS surveys. The six galaxies are visually selected to have clumpy substructures with distinct colors over rest frame 3600–4100 Å, including a red, dominant stellar core that is close to their stellar-light centroids. With 23-filter photometry from the Hubble Space Telescope to JWST, we measure the stellar-population properties of individual structural components via spectral energy distribution fitting using PROSPECTOR. We find that the central stellar cores are $\gtrsim 2$ times more massive than the Toomre mass, indicating they may not form via single in situ fragmentation. The stellar cores have stellar ages of 0.4–0.7 Gyr that are similar to the timescale of clump inward migration due to dynamical friction, suggesting that they likely instead formed through the coalescence of giant stellar clumps. While they have not yet quenched, the six galaxies are below the star-forming main sequence by 0.2–0.7 dex. Within each galaxy, we find that the specific star formation rate is lower in the central stellar core, and the stellar-mass surface density of the core is already similar to quenched galaxies of the same masses and redshifts. Meanwhile, the stellar ages of the cores are either comparable to or younger than the extended, smooth parts of the galaxies. Our findings are consistent with model predictions of the gas-rich compaction scenario for the buildup of galaxies’ central regions at high redshifts. We are likely witnessing the coeval formation of dense central cores, along with the onset of galaxy-wide quenching at $z > 3$.

*Unified Astronomy Thesaurus concepts: Galaxy formation (595); Galaxy evolution (594); Galaxy quenching (2040); Galaxy structure (622); High-redshift galaxies (734)*

## 1. Introduction

Over the last two decades, deep imaging surveys with the Hubble Space Telescope (HST), such as GOODS (M. Giavalisco et al. 2004), CANDELS (A. M. Koekemoer et al. 2011; N. A. Grogin et al. 2011), and Hubble Frontier Fields

(J. M. Lotz et al. 2017), have obtained more than 100,000 rest-frame optical images for mass complete samples of galaxies up to redshift $z=3$. These observations have not only greatly advanced our understanding of galaxies’ structural evolution across cosmic time, but also raised a number of new questions.

It immediately becomes apparent after analyzing those HST images of $1 < z < 3$ galaxies that they have distinct structural and morphological properties. Relative to galaxies in the local Universe, massive galaxies at $z \sim 2$ are on average smaller and

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA; zhiyuanji@arizona.edu</footnote>

<footnote n="2">²NSF's National Optical-Infrared Astronomy Research Laboratory, 950 North Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="3">³Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK</footnote>

<footnote n="4">⁴Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK</footnote>

<footnote n="5">⁵Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="6">⁶Kavli Institute for Particle Astrophysics and Cosmology and Department of Physics, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="7">⁷Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="8">⁸Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge MA 02138, USA</footnote>

<footnote n="9">⁹Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter Street, Madison, WI 53706, USA</footnote>

<footnote n="10">¹⁰Department of Physics and Astronomy, The Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA</footnote>

<footnote n="11">¹¹Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

<footnote n="12">¹²Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="13">¹³Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield AL10 9AB, UK</footnote>

<footnote n="14">¹⁴Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="15">¹⁵NRC Herzberg, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada</footnote>

<footnote n="16">¹⁶Department for Astrophysical and Planetary Science, University of Colorado, Boulder, CO 80309, USA</footnote>

<footnote n="17">¹⁷European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain</footnote>

<footnote n="18">¹⁸European Space Agency, ESA/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands</footnote>

<footnote n="19">¹⁹School of Physics, University of Melbourne, Parkville 3010, VIC, Australia</footnote>

<footnote n="20">²⁰ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia</footnote>

<footnote>Received 2023 May 29; revised 2024 July 24; accepted 2024 August 7; published 2024 October 9</footnote>
denser (e.g., H. C. Ferguson et al. 2004; A. van der Wel et al. 2014; T. Shibuya et al. 2015; G. Barro et al. 2017), and their morphologies are irregular and often show large-scale asymmetries (e.g., L. L. Cowie et al. 1995; S. van den Bergh et al. 1996). This latter morphological feature is also present in $z \sim 0$ galaxies, but it has been predominately observed in galaxy mergers (e.g., C. J. Conselice et al. 2000; J. M. Lotz et al. 2004). Despite the fact that galaxy major merger rate increases with redshift (C. J. Conselice 2014), spatially resolved spectroscopy revealed that at least half of the $z \sim 2$ star-forming galaxies have a rotationally supported disk, arguing against them being interacting systems (e.g., B. J. Weiner et al. 2006; N. M. Förster Schreiber et al. 2006, 2009, 2018; R. Genzel et al. 2006; K. L. Shapiro et al. 2008; D. R. Law et al. 2009; S. A. Kassin et al. 2012; R. C. Simons et al. 2017; E. Wisnioski et al. 2019; N. M. Förster Schreiber & S. Wuyts 2020). Moreover, most of the $z \sim 2$ star-forming disks show several giant, kiloparsec-scale stellar clumps (e.g., B. G. Elmegreen & D. M. Elmegreen 2005; D. M. Elmegreen et al. 2007, 2009; R. Genzel et al. 2008, 2011, 2014; N. M. Förster Schreiber et al. 2011; Y. Guo et al. 2012, 2015), a morphological characteristic rarely seen in the local Universe. Understanding this substantial change in galaxies' structures across cosmic time is the key to constraining models not only of galaxy evolution but also of the interplay between baryonic and dark matter (e.g., R. Genzel et al. 2020).

Another important process tied to structural transformation is galaxy quenching—the physical processes that control the transition from widespread star formation to quiescence in galaxies. Since $z=3$, star-forming galaxies showed a gradual transformation in their structures from disturbed to normal, smooth disks (M. Huertas-Company et al. 2015), while the fraction of clumpy star-forming galaxies decreased from $\approx 50\%$ to $\approx 5\%$ (e.g., K. L. Murata et al. 2014; Y. Guo et al. 2015). In the meantime, there was a rapid buildup of quiescent galaxies (e.g., O. Ilbert et al. 2010; A. Muzzin et al. 2013; O. Ilbert et al. 2013; A. R. Tomczak et al. 2014), which notably have much larger stellar-mass surface densities than star-forming galaxies of the same masses and redshifts (e.g., G. Kauffmann et al. 2003; J. Brinchmann et al. 2004; M. Franx et al. 2008; E. Cheung et al. 2012; M. Mosleh et al. 2017). Is there a causal link between structural transformations and quenching? While a physical link between the two phenomena has been seen in cosmological simulations (A. Zolotov et al. 2015; S. Wellons et al. 2015; S. Tacchella et al. 2016), this question still eludes us observationally, because the structure of galaxies is only known at the time of observation, which makes empirical constraints on the relative timing sequence of the two events difficult (e.g., K. Bundy et al. 2010; A. van der Wel et al. 2011; V. A. Bruce et al. 2012; A. B. Newman et al. 2015, 2018; S. Toft et al. 2017; Z. Ji & M. Giavalisco 2022, 2023).

There is growing evidence that, as they head toward quiescence, galaxies develop a dense stellar core, possibly a bulge (e.g., G. Kauffmann et al. 2003; J. Brinchmann et al. 2004; E. Cheung et al. 2012; C. C. Williams et al. 2014; E. Nelson et al. 2014; P. Lang et al. 2014; V. A. Bruce et al. 2014; S. Tacchella et al. 2015; G. Barro et al. 2017; K. E. Whitaker et al. 2017; K. A. Suess et al. 2021; P. Dimauro et al. 2022; Z. Ji & M. Giavalisco 2023), but the exact physical mechanism responsible for this buildup of central regions remains unknown. While major mergers and weak,

nonaxisymmetric instabilities (e.g., bars, spiral arms) are believed to drive the development of central structures in $z \sim 0$ disks (J. Kormendy & R.C.J. Kennicutt 2004), the physical conditions, hydrodynamical instabilities in particular, are very different in high-redshift galaxies. One mechanism to effectively and efficiently grow the central regions of a galaxy at high redshift is actually a class of processes, generically referred to as “wet/gas-rich compaction” whose main feature is highly dissipative gas accretion toward the galactic center, which, in turn, promotes enhanced activity of star formation at a higher rate relative to the average star formation rate (SFR; A. Dekel et al. 2009). Consequently, right after a compaction event, a galaxy will have a younger center relative to the other smooth parts of the galaxy. In the absence of further gas inflow, this leads to the central gas depletion and finally the quenching of central star formation (e.g., S. Tacchella et al. 2016). Hydrodynamical simulations with adequate resolution to model relevant gaseous physics show that gas-rich major mergers (A. Zolotov et al. 2015; S. Inoue et al. 2016), collisions of counterrotating gas streams (M. Danovich et al. 2015), and violent disk instability (VDI) are all associated with gas-rich compaction.

Regarding VDI in particular, it can be understood in a framework of Toomre instability (A. Toomre 1964; J. Binney & S. Tremaine 2008), where a rotating gaseous disk becomes unstable when$^{21}$

$$Q = \frac{\sigma\kappa}{\pi G \Sigma} < Q_c \sim 1 \quad (1)$$

where $\sigma$ is local gas velocity dispersion, $\kappa$ is epicyclic frequency, and $\Sigma$ is local gas surface density. Observations found that star-forming galaxies at cosmic noon are characterized by large $\sigma$, with a typical level of $50 \text{ km s}^{-1}$ (e.g., R. Genzel et al. 2008; N. M. Förster Schreiber et al. 2009; R. Genzel et al. 2011) compared to only a few kilometers per second observed for nearby spiral galaxies. In addition, high-redshift galaxies have higher gas fractions, with a typical value of $\approx 50\%$ compared to a few percent to 10% at $z \sim 0$ (L. J. Tacconi et al. 2020). With these properties, high-redshift disks have $Q \sim 1$ (i.e., the entire disk is unstable/marginally stable), while low-redshift disks are largely stable with $Q > 1$. These lead to two direct consequences in high-redshift galaxies. First, high-redshift disks are subject to VDI that lead to the formation of giant stellar clumps. Second, the timescale for the inward migration of clumps toward the gravitational center$^{22}$ is much shorter at high redshifts than in the nearby Universe (A. Dekel et al. 2009; R. Genzel et al. 2011).

A popular hypothesis is thus that massive stellar cores in high-redshift galaxies are the result of unstable, gas-rich disks, which first fragment into giant stellar clumps that then rapidly migrate inward and finally coalesce to form dense stellar cores (e.g., M. Noguchi 1999; A. Immeli et al. 2004; F. Bournaud et al. 2007, 2014; B. G. Elmegreen et al. 2008; A. Dekel et al. 2009; D. Ceverino et al. 2010). On the one hand, the observed radial color gradients of star-forming clumps seem to support this hypothesis (N. M. Förster Schreiber et al. 2011; Y. Guo et al. 2012, 2015; T. Shibuya et al. 2016; E. Soto et al. 2017;

<footnote>$^{21}$ For thick disks, $Q_c$ is about 0.7 rather than 1 (M. Behrendt et al. 2015).</footnote>

<footnote>$^{22}$ The timescale is set by dynamical friction; hence, it is proportional to the square of the ratio of disk rotation velocity divided by velocity dispersion, i.e., $(v_{rot}/\sigma)^2$.</footnote>
<table><caption>**Table 1**

Physical Properties of the Final Sample Presented in This Work</caption><thead><tr><th>ID</th><th>R.A.<br>(Deg)</th><th>Decl.<br>(Deg)</th><th>Redshift</th><th>Redshift Ref.<sup>a</sup></th><th>log M<sub>*</sub><sup>b</sup><br>(M<sub>☉</sub>)</th><th>SFR<sup>b</sup><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>Stellar Age<sup>b,c</sup><br>(Gyr)</th></tr></thead><tbody><tr><td>JADES-JEMS-13396</td><td>53.14921</td><td>-27.79158</td><td>3.559</td><td>MUSE</td><td>10.00<sup>+0.02</sup><sub>-0.04</sub></td><td>13.5<sup>+0.8</sup><sub>-0.6</sub></td><td>0.8<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>JADES-JEMS-15157</td><td>53.13943</td><td>-27.78009</td><td>3.591</td><td>MUSE and FRESCO<sup>d</sup></td><td>9.36<sup>+0.09</sup><sub>-0.11</sub></td><td>1.3<sup>+1.0</sup><sub>-0.4</sub></td><td>0.4<sup>+0.2</sup><sub>-0.1</sub></td></tr><tr><td>JADES-JEMS-6885</td><td>53.12832</td><td>-27.84573</td><td>3.698</td><td>FRESCO</td><td>9.87<sup>+0.03</sup><sub>-0.07</sub></td><td>7.4<sup>+1.1</sup><sub>-0.9</sub></td><td>0.4<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>JADES-JEMS-16296</td><td>53.15651</td><td>-27.77227</td><td>3.53</td><td>Photometric</td><td>9.65<sup>+0.05</sup><sub>-0.06</sub></td><td>5.4<sup>+1.2</sup><sub>-0.9</sub></td><td>0.5<sup>+0.3</sup><sub>-0.2</sub></td></tr><tr><td>JADES-JEMS-14436</td><td>53.13306</td><td>-27.78456</td><td>3.67</td><td>Photometric</td><td>10.34<sup>+0.07</sup><sub>-0.09</sub></td><td>11.8<sup>+2.5</sup><sub>-1.9</sub></td><td>0.8<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>JADES-JEMS-11059</td><td>53.13484</td><td>-27.80887</td><td>4.07</td><td>Photometric</td><td>9.59<sup>+0.08</sup><sub>-0.07</sub></td><td>9.6<sup>+0.7</sup><sub>-0.6</sub></td><td>0.4<sup>+0.1</sup><sub>-0.2</sub></td></tr></tbody></table>

**Notes.**

<sup>a</sup> Redshifts are from either the MUSE/HUDF Data Release 2 (R. Bacon et al. 2022), or [S III]λ9531 emission line detected in the NIRCam WFSS spectra (F. Sun 2024, private communication) obtained by the FRESCO survey (P. A. Oesch et al. 2023), or photometric redshifts when spectroscopic ones are not available.

<sup>b</sup> Reported values are from the PROSPECTOR spectral energy distribution (SED) fitting with our fiducial model (Section 2.4).

<sup>c</sup> Mass-weighted stellar age.

<sup>d</sup> The redshift of this galaxy derived using the [S III]λ9531 in its FRESCO spectrum is slightly different from using the Lyα in its MUSE spectrum, with the former being 3.591 and latter being 3.605.

N. Mandelker et al. 2017). On the other, studies of gas kinematics find evidence that strong outflows can launch from the clump sites (R. Genzel et al. 2011), arguing that clumps might be self-disrupted by strong stellar feedback, which can significantly reduce the growth rate of central stellar cores (N. Murray et al. 2010; P. F. Hopkins et al. 2012; S. Genel et al. 2012; T. Buck et al. 2017). Thus, it remains unknown to what extent VDIs drive the growth of central regions in high-redshift galaxies.

The James Webb Space Telescope (JWST; J. P. Gardner et al. 2023) is revolutionary for studying the detailed physics of structural transformations and quenching at cosmic noon. For the first time, we can get images at rest-frame optical and near-IR (NIR) wavelengths, with a sensitivity that is not attainable with HST, for galaxies beyond $z > 3$, closer to the time of early core/bulge formation and quenching. This guarantees more direct constraints on the underlying physics. With subkiloparsec spatial resolution, JWST’s Near Infrared Camera (NIRCam; M. J. Rieke et al. 2023b) can resolve a physical scale comparable to the typical Toomre length of gaseous disks at cosmic noon (e.g., A. Escala & R. B. Larson 2008; A. Dekel et al. 2009; D. M. Elmegreen et al. 2009; R. Genzel et al. 2011), promising to identify individual stellar substructures. Together with the ancillary HST data, imaging that covers truly panchromatic swathes from rest-frame UV through NIR enable us to put unprecedented constraints on the detailed, spatially resolved mass assembly history of galaxies at high redshifts.

In this paper, we present the spatially resolved stellar populations in six galaxies at $3 < z < 4.5$. This redshift range represents one of the key epochs of the Universe when quenching becomes a significant process (A. Muzzin et al. 2013; C. Schreiber et al. 2018; E. Merlin et al. 2019; A. Shahidi et al. 2020; A. C. Carnall et al. 2023b, 2023a), which can provide unique insight into the role that substructures among massive galaxies may play in the transition from star-forming disks to bulge-dominated quiescent objects. The galaxies are visually selected to have bright and red central regions. We show that these galaxies are very likely in the act of developing their stellar cores/bulges. We discuss their implications on structural transformations and quenching at cosmic noon. Throughout this paper, we adopt the Λ cold dark matter cosmology with Planck Collaboration et al. (2020) parameters, i.e., $\Omega_m = 0.315$ and $h = H_0/(100 \text{ km s}^{-1} \text{ Mpc}^{-1}) = 0.673$.

## 2. Observations and Analysis

### 2.1. Sample Selection

This work is a pilot study focused only on six galaxies at $3 < z < 4.5$ having visually identified prominent, red stellar cores at rest-frame 4000 Å. Basic information of the six galaxies is present in Table 1. Utilizing the power of multiband JWST+HST images, we provide a spatially resolved view of the stellar-population properties of these six galaxies, aiming at shedding light on the physical mechanisms responsible for the formation of central spheroidal structure in massive galaxies at cosmic noon.

The six galaxies are drawn from an initial sample of HST/F160W-selected galaxies with $3 < z < 4.5$ from the latest CANDELS/GOODS-S catalog obtained by the ASTRODEEP project (E. Merlin et al. 2021). By selection, we first retain 361 galaxies in the footprint of the JWST Extragalactic Medium-band Survey (JEMS; ≈10 arcmin²; C. C. Williams et al. 2023b, JWST Cycle 1, PID: 1963), which also appear in the footprint of the JWST Advanced Deep Extragalactic Survey (JADES; D. J. Eisenstein et al. 2017), a collaborative program between the NIRCam and NIRSpec GTO teams. These yield a total of 14-filter NIRCam imaging at $\lambda \approx 0.8-5 \mu$m (Section 2.2). Together with the ancillary data from HST, this high angular resolution imaging data across a total of 23 filters enables us to perform detailed stellar-population analysis in these galaxies. With the new photometry from NIRCam imaging (Section 2.2), we then run spectral energy distribution (SED) fitting with the code PROSPECTOR (B. D. Johnson et al. 2021) using the fiducial model described in Section 2.4, during which we set redshift as a free parameter. Finally, we only retain 301 galaxies as the parent sample whose best-fit redshifts from our new SED fitting are still in the range of $3 < z_{\text{phot}} < 4.5$.

For each one of the galaxies in the parent sample, we visually check the color maps produced by its NIRCam F150W, F182M, and F210M images, which probe rest-frame light at ~3000, 4000, and 4500 Å, respectively. This wavelength range, i.e., around the 4000 Å break, is chosen because it contains key information about stellar-population properties, especially stellar ages (e.g., A.G. Bruzual 1983; G. Kauffmann et al. 2003). Our visual identification finds 37 (out of 301) galaxies having substructures with distinct F150W/F182M/F210M colors. As Figure 1 shows, these 37
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The parent sample of $3 < z < 4.5$ galaxies from JEMS. The x-axis shows the stellar mass from our PROSPECTOR SED fitting (Section 2.4), the y-axis shows the half-light radius derived directly from the F210M surface brightness profile (no point-spread function (PSF)-convolved analysis included). Red circles show the 37 galaxies with visually identified substructures having distinct rest-frame colors at 4000 Å. The golden hexagons show the final sample of six galaxies with red, massive cores presented in this work (see Section 2.1 for the detailed sample selection).</figcaption></figure>

galaxies are at the massive end of our JEMS parent sample, i.e., the majority of them have stellar masses >$10^9 M_\odot$, which is expected because the visual selection is biased toward bright substructures in relatively bright galaxies. Among 36 galaxies in the parent sample with masses >$10^{9.5} M_\odot$, 22 of them are visually identified to have substructures, corresponding to a fraction of 61% ± 17%. Interestingly, at face value this fraction is in quantitative agreement with earlier HST studies of the fraction of galaxies with UV clumps at slightly lower redshifts $2.5 < z < 3$ (e.g., T. Shibuya et al. 2016; Z. Sattari et al. 2023). While we caution that the methods used to identify clumps can significantly affect the results, we defer more detailed investigation of this impact to a future work with larger samples and more uniform selection methods.

As Figure 2 shows, our analysis finds that the 37 galaxies with visually identified substructures have diverse spatially resolved stellar populations, representing different stages of galaxy structural transformations. This motivates the present paper and upcoming papers in this series, where we address relevant physical mechanisms associated with those different stages of structural transformations.

In this first paper, we only focus on a group of six galaxies whose main characteristic is to have a bright and red (at rest-frame ~4000 Å) stellar core, which we refer to as “C1” in the following. C1 is close to the center of stellar-light distribution, and it is often accompanied by additional off-center, minor clumpy substructures (named “C2” and “C3,” if present). Figures 3 and 4 show their 23-filter cutouts. Our analysis below shows that C1 is more massive than what one would expect from a single collapse due to the Toomre instability, and has a star formation history (SFH) that shows a recent major star formation episode started ~0.5–1 Gyr ago that peaked at ~0.1–0.3 Gyr prior to the time of observation. The formation of the stellar core C1 and its impact, if any, on star formation will be the main focus of this paper.

Finally, we note that in the future, instead of using visual selection, we will use a more uniform selection method to study

resolved properties of the more general population of galaxies at $3 < z < 4.5$ with deep JADES NIRCam imaging.

## 2.2. Observations and Data Reduction

The six galaxies have been observed with deep$^{23}$ imaging in 14 JWST/NIRCam filters. These include nine-filter images (F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, and F444W) acquired as part in the GOODS-S portion of JADES, and five additional medium-filter images (F182M, F210M, F430M, F460M, and F480M) acquired by JEMS.

The NIRCam observations of all 14 filters are reduced using a consistent method, which has been briefly summarized in B. E. Robertson et al. (2022) and will be presented in detail in M. Rieke et al. (2023a) and S. Tacchella et al. (2023, in preparation). In short, we process the raw images with the JWST Calibration Pipeline v1.8.1 (H. Bushouse et al. 2022) with the CRDS pipeline mapping (pmap) context 1009. To perform the detector-level corrections and obtain count-rate images, we run the Stage 1 of the JWST pipeline with default parameters, during which we also mask and correct for the “snowballs” effect caused by charge deposition arising from cosmic ray hits. We then run the Stage 2 of the JWST pipeline, with default parameters, to perform the flat-fielding and flux calibration, after which we perform several custom corrections including the removal of $1/f$ noise and the subtraction of 2D background from the images, and the correction for the “wisp” features in the short-wavelength (SW, i.e., filters bluer than F277W) images. Afterwards, for individual exposures of a given visit, we match sources to the World Coordinate System (WCS) of a reference catalog constructed from HST F160W mosaics in the GOODS-S field with astrometry tied to Gaia-EDR3 (Gaia Collaboration et al. 2021), calculate astrometric corrections, and then run the Stage 3 of the JWST pipeline to combine together all of those individual exposures. Finally, we combine individual visit-level mosaics together to create the final mosaic, where we choose a pixel scale of 0.03 pixel$^{-1}$ for all filters.

The ancillary HST data used in this work are taken from the latest Hubble Legacy Fields Data Release$^{24}$ in the GOODS-S field (G. Illingworth et al. 2016; K. E. Whitaker et al. 2019), including the imaging data of the ACS/WFC F435W, F606W, F775W, F850LP, and F814W filters, and of the WFC3 F105W, F125W, F140W, and F160W filters. For all of those HST images, we have (1) tied them to the same astrometry and (2) resampled them to the same pixel scale as those of our NIRCam images (B. Robertson et al. 2023, in preparation).

## 2.3. Photometry

For each galaxy, we perform aperture-matched photometry for (1) integrated flux, (2) individual substructures, and (3) a smooth component. The integrated flux is defined as total flux in the pixel aperture generated from a galaxy’s segmentation map generated by JADES (B. E. Robertson et al. 2022). For fluxes of individual substructures, we use customized elliptical apertures that are illustrated in the RGB images in Figures 6 and 7. Finally, the smooth component is defined as integrated flux minus coadded fluxes of individual substructures.

<footnote>$^{23}$ The sensitivity for a 5σ detection using a $r = 0''15$ circular aperture is ~30 mag in the F200W filter (B. E. Robertson et al. 2022).</footnote>

<footnote>$^{24}$ https://archive.stsci.edu/prepds/hlf/</footnote>
<figure><img src="image_3.png" /></figure>

**Figure 2.** Examples of galaxies with different types of substructures identified by our visual selection (see Section 2.1). The RGB cutouts shown here are made with JWST/NIRCam F150W (B), F182M (G), and F210M (R) images. The pink circle marks the angular resolution (FWHM) of NIRCam/F210M imaging. Our visual identifications find a sample of 37 galaxies having substructures with distinct colors at rest-frame 4000 Å. These galaxies show diverse substructures, and can be generally categorized into three groups. The left panel shows an example galaxy with a red, prominent stellar core. The middle panel shows an example galaxy with a dusty core, where in the RGB (made with images at rest-4000 Å) cutout, a flux deficit is seen near the center of the NIRCam/F480M image (rest-1 µm, starlight). The right panel shows an example galaxy with a compact blue center. In this paper, we will only focus on six galaxies having visually identified prominent, red stellar cores at rest-frame 4000 Å (similar to the example galaxy shown in the left panel), while the latter (middle and right) two types of galaxies will be studied in forthcoming papers.

Point-spread functions (PSFs) of all images are first homogenized to that of NIRCam/F444W. To do so, in the footprint of JEMS, we visually identify 13 isolated, bright stars that are unsaturated across all filters from HST through JWST. We then build effective PSFs (ePSFs) following the method from J. Anderson & I. R. King (2000). For NIRCam filters, we have also generated model PSFs (mPSFs) using the package WEBBPSF (v1.1.1, M. D. Perrin et al. 2012, 2014), during which we also took into account the observational (e.g., source location on detectors, optical path difference) and data reduction (e.g., mosaicking) effects. As shown in Appendix A, we find excellent agreement between ePSFs and mPSFs, with a typical difference of ≲1%. We then homogenize ePSFs using the package PYPHER, which calculates the convolution kernel based on Wiener filtering with a regularization factor (A. Boucaud et al. 2016). Finally, we perform aperture-matched photometry in the PSF-matched images.

To estimate photometric uncertainties, we follow I. Labbé et al. (2005) to first use randomly placed apertures in empty (i.e., no sources) regions in the images to estimate the growth of noise as a function of the linear size of aperture. We then fit the function using the parameterization of R. Quadri et al. (2007) and finally estimate photometric uncertainties using Equation 5 in K. E. Whitaker et al. (2011). This procedure will be presented in detail in B. Robertson et al. (2023, in preparation).

## 2.4. Measuring Stellar-population Properties

We model panchromatic SEDs to derive stellar-population properties of individual structural components, with the emphasis on robustly measuring their stellar masses, ages, and SFHs. For each one of the galaxies, we perform two rounds of SED fitting. The first round aims at checking the physical association of all structural components, as opposed to the case when sources at different distances are projected onto the regions of sky in close proximity of each other. During the SED fitting of this round, we set redshift as a free parameter with a flat prior of $z \in (0, 10)$. We then check posterior distributions of redshift, and consider structural components as being physically linked when the difference in the median of their redshift posteriors is ≲0.1. Afterwards, we fix redshifts of all components to either the spectroscopic redshift if available

or the best-fit value that we derive from the first-round SED fitting using the galaxy's integrated photometry,<fnref n="25" /> and run the second-round SED fitting. The fitting results are shown in Figures 6 and 7, and will be discussed in Section 3.1. Below, we detail the model assumptions of our SED fitting.

### 2.4.1. Basic Setup of the Prospector Fitting

The SED fitting code PROSPECTOR (B. D. Johnson et al. 2021) is elected to use in this work. It is built upon a fully Bayesian framework that makes it possible to fit galaxies' SEDs with complex models of stellar population synthesis.

Regarding the basic setup of PROSPECTOR fitting, we adopt the Flexible Stellar Population Synthesis (FSPS) code (C. Conroy et al. 2009; C. Conroy & J. E. Gunn 2010) where we use the stellar isochrone libraries MIST (A. Dotter 2016; J. Choi et al. 2016) and the stellar spectral libraries MILES (J. Falcón-Barroso et al. 2011). During the modeling, we use the Markov Chain Monte Carlo (MCMC) sampling code DYNESTY (J. S. Speagle 2020), which adopts the nested sampling procedure (J. Skilling 2004). We assume the P. Kroupa (2001) initial mass function (IMF) for consistency with the IMF used in the N. Byler et al. (2017) nebular continuum and line emission model, which we adopt in our SED modeling. We adopt the P. Madau (1995) intergalactic medium absorption model.

We set the stellar metallicity as a free parameter and assume a flat prior in logarithmic space $log(Z*/Z_\odot) \in (-2, 0.19)$, where $Z_\odot = 0.0142$ is solar metallicity. The upper limit of the prior is chosen because it is the highest metallicity that the MILES library has. We also leave the gas-phase metallicity $Z_{gas}$ and ionization parameter $U$ as free parameters during the fitting, where we use flat priors of $log(Z_{gas}/Z_\odot) \in (-2, 0.5)$ and $log U \in (-4, 1)$.

Following S. Tacchella et al. (2022a), we assume a two-component dust attenuation model where the dust attenuation of nebular emission and young stellar populations, and of old stellar populations, are treated differently (S. Charlot & S. M. Fall 2000). For stellar populations older than 10 Myr,

<footnote n="25">²⁵ We notice that, for galaxies with spectroscopic redshifts, our photometric redshifts are highly consistent with the spectroscopic ones, with a typical difference of $|z_{spec} - z_{phot}| < 0.1$, thanks to the high-S/N photometric coverage around both the Lyman Break and Balmer Break.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Panchromatic HST and JWST/NIRCam images of the three galaxies with spectroscopically confirmed redshifts $z_{spec}$. The RGB image shown on the left is made with JWST/NIRCam F150W (B), F182M (G), and F210M (R) images. We use a linear color scale for the RGB cutout to better show the substructures. The pink circle at the bottom-left marks the FWHM of the NIRCam/F210M PSF (Appendix A). Small panels show cutouts of 23 individual filters covering observed $\lambda \approx 0.35 -5\mu m$ included in the analysis of this work, where the newly obtained 14-filter NIRCam images are labeled with a hexagon symbol at the bottom right. The last small panel shows the segmentation map. The cutouts of individual filters are shown in a log color scale, highlighting the smooth component of the galaxies.</figcaption></figure>

we assume the dust attenuation using the parameterization from S. Noll et al. (2009), i.e.,

$$ \tau_{\text{dust},2}(\lambda) = \frac{\hat{\tau}_{\text{dust},2}}{4.05} (k'(\lambda) + D(\lambda)) \left( \frac{\lambda}{5500\AA} \right)^{\delta} \quad (2) $$

where $\hat{\tau}_{\text{dust},2}$ and $\delta$ correspond to dust2 and dust_index in the FSPS, $k'(\lambda)$ is the D. Calzetti et al. (2000) dust attenuation

law, and $D(\lambda)$ is the Lorentzian-like profile used to describe the UV dust bump at 2175 Å. Instead of setting $D(\lambda)$ as a free parameter, we tie it to $\delta$ following the results of M. Kriek & C. Conroy (2013). Therefore, the dust attenuation for the older stellar populations only has two free parameters, $\delta$ and $\hat{\tau}_{\text{dust},2}$. We assume flat priors for both parameters, i.e., $\delta \in (-1, 0.4)$ and $\hat{\tau}_{\text{dust},2} \in (0, 10)$. For stellar populations younger than
<figure><img src="image_5.png" /><figcaption>Figure 4. Similar to Figure 3, but for the remaining three galaxies with photometric redshifts $z_{phot}$.</figcaption></figure>

10 Myr, we assume the same dust attenuation law as for the nebular emission, which has a functional form of

$$ \tau_{\text{dust},1}(\lambda) = \hat{\tau}_{\text{dust},1} \left( \frac{\lambda}{5500\AA} \right)^{-1}, \quad (3) $$

where $\hat{\tau}_{\text{dust},1}$ corresponds to the parameter dust1 in the FSPS. Instead of modeling $\hat{\tau}_{\text{dust},1}$ as an independent free parameter, we tie it to $\hat{\tau}_{\text{dust},2}$ and model their ratio using a clipped normal prior centered at 1, with a width of 0.3 and in the range of $\hat{\tau}_{\text{dust},1}/\hat{\tau}_{\text{dust},2} \in (0, 2)$.

We also include active galactic nucleus (AGN) dust torus templates from M. Nenkova et al. (2008a) and M. Nenkova et al. (2008b), but we stress that the purpose of including this AGN component is not to quantitatively constrain the AGN strength, since the data coverage we currently have does not have the power to do so. Instead, the purpose simply is to marginalize over it in order to check how the AGN component can affect the fitting results. This adds two free parameters: the ratio of bolometric luminosity from the galaxy divided by that from the AGN ($f_{AGN}$), and the optical depth of clumps in AGN dust torus at 5500 Å ($\tau_{AGN}$). We assume flat priors in logarithmic space for both parameters, i.e., $f_{AGN} \in (10^{-5}, 3)$
and $\tau_{AGN} \in (5, 150)$. We find that $f_{AGN}$ is consistent with zero for all sample galaxies, which is in line with the fact that these galaxies are not classified as AGNs using other different selection methods, including X-ray (B. Luo et al. 2017), IRAC colors (Z. Ji et al. 2022), and multiwavelength selections (J. Lyu et al. 2022).

### 2.4.2. Reconstructing Star Formation Histories

One key feature of PROSPECTOR is that it allows flexible parameterizations of galaxies' SFHs, both in parametric and nonparametric forms. This has been demonstrated as crucial to reducing systematic biases in the measurements of stellar mass, SFR, and stellar age (J. Leja et al. 2019). Several recent studies have further shown that statistically, PROSPECTOR is capable of reconstructing high-fidelity nonparametric SFHs for synthetic galaxies generated from cosmological simulations (J. Leja et al. 2019; B. D. Johnson et al. 2021; S. Tacchella et al. 2022a; Z. Ji & M. Giavalisco 2022).

For the fiducial SFH, we assume a piecewise, nonparametric form composed of $N_{SFH} = 9$ lookback time bins ($t_{lookback}$, i.e., the time prior to the time of observation), where SFR is constant in each bin. Among the nine lookback time bins, the first two bins are fixed to be 0–30 and 30–100 Myr in order to capture the recent star formation activity; the last bin is assumed to be 0.85$t_H$–$t_H$ where $t_H$ is the Hubble Time at the time of observation; and the remaining six bins are evenly spaced in logarithmic space between 100 Myr and 0.85$t_H$. As has been extensively tested by J. Leja et al. (2019) using mock observations of simulated galaxies (see their Figure 15), the recovered physical properties are largely insensitive to $N_{SFH}$ when it is greater than 5. This is also found later through our analysis presented in Section 3.1 and Figure 5.

The procedure above to reconstruct nonparametric SFHs potentially suffers from overfitting problems, because it includes “more bins than the data warrant” (J. Leja et al. 2019). An effective solution to mitigate the issue is to choose a prior weight for physically plausible forms, as opposed to letting the SFHs have a fully arbitrary shape (A. C. Carnall et al. 2019a; J. Leja et al. 2019). Our fiducial measures use the continuity prior that has been demonstrated to work well across various galaxy types (J. Leja et al. 2019). Recent studies have noted the systematics introduced by assumed priors to reconstructed SFHs (K. A. Suess et al. 2022; S. Tacchella et al. 2022b; Z. Ji & M. Giavalisco 2022). To check possible systematics in our measurements, we therefore also experiment using another well-tested Dirichlet prior (J. Leja et al. 2017). In addition, we have also tested our results using the parametric delayed-tau model. As we show in Appendix C, the results of this work remain qualitatively unchanged using these different methods of SFH reconstructions. In the remainder of main text, our discussion thus will only focus on measurements from the continuity prior.

## 2.5. Measuring Morphological Properties of the Substructures

To estimate intrinsic stellar masses (Section 3.2) and stellar-mass surface densities (Section 4.2) of the clumpy substructures, we perform PSF-convolved morphological fitting. For each one of the galaxies, we simultaneously model the 2D light distributions of all clumpy substructures with the fitting code FORCEPHO (B. Johnson et al. 2023, in preparation), which adopts a fully MCMC framework and performs the fitting on

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparisons of the stellar-mass (top panel) and star formation rate (SFR; bottom panel) measures. The x-axis shows the measure using integrated flux, and the y-axis shows the results by coadding the measures of individual structural components (i.e., substructures+smooth) together. We plot all 37 galaxies showing distinct substructures around the rest-frame 4000 Å (Section 2.1), and the final sample of six galaxies presented in this paper are marked with the red star symbols. We see excellent agreement between the two measurements.</figcaption></figure>

the basis of individual NIRCam exposures, as opposed to performing the fitting in a drizzled image (i.e., a mosaic). In this way, parameter uncertainties are better estimated because correlations of adjacent image pixels, and, more importantly here given the complex morphologies of the galaxies, correlations of fluxes from neighboring substructures are naturally taken into account (see B. E. Robertson et al. 2022 and S. Tacchella et al. 2023 for a brief discussion).

A redder NIRCam filter generally better probes the stellar-mass distribution. Meanwhile, a better angular resolution also helps to better separate fluxes from individual substructures. The combination of the poorer imaging angular resolution of longer-wavelength filters and the younger stellar populations of off-center clumps (i.e., their SEDs are bluer) makes it increasingly difficult to morphologically decompose individual substructures using images of longer-wavelength filters. We thus perform the FORCEPHO analysis using the reddest NIRCam's SW filter we have, i.e., F210M, which probes rest-frame ~5000 Å of the six galaxies. Finally, we model each clumpy substructure using a single Sérsic profile.

Figure 8 shows results of the FORCEPHO fitting. As the residual map shows, each one of the substructures is well fit by a single Sérsic profile. The χ map also confirms the presence of
<table><caption>**Table 2**
Properties of Spatially Resolved Stellar Populations</caption><thead><tr><th>ID</th><th>Redshift</th><th>Component</th><th>F<sub>v</sub><sup>F150Wa</sup><br>(nJy)</th><th>log M<sub>*</sub><sup>b</sup><br>(M<sub>☉</sub>)</th><th>SFR<sup>b,c</sup><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>Stellar Age<sup>b</sup><br>(Gyr)</th><th>sSFR<sup>b</sup><br>(Gyr<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="3">JADES-JEMS-13396</td><td rowspan="3">3.559</td><td>C1</td><td>24.6 ± 0.8</td><td>9.40 <sup>+0.03</sup><sub>-0.04</sub></td><td>2.7 <sup>+0.2</sup><sub>-0.3</sub></td><td>0.7 <sup>+0.1</sup><sub>-0.2</sub></td><td>1.07 <sup>+0.10</sup><sub>-0.15</sub></td></tr><tr><td>Smooth</td><td>212.8 ± 11.4</td><td>9.88 <sup>+0.06</sup><sub>-0.13</sub></td><td>13.4 <sup>+2.3</sup><sub>-1.6</sub></td><td>0.6 <sup>+0.1</sup><sub>-0.2</sub></td><td>1.76 <sup>+0.38</sup><sub>-0.56</sub></td></tr><tr><td>Integrated</td><td>237.4 ± 11.7</td><td>10.00 <sup>+0.02</sup><sub>-0.04</sub></td><td>13.5 <sup>+0.8</sup><sub>-0.6</sub></td><td>0.8 <sup>+0.1</sup><sub>-0.1</sub></td><td>1.35 <sup>+0.10</sup><sub>-0.13</sub></td></tr><tr><td rowspan="6">JADES-JEMS-15157</td><td rowspan="6">3.591</td><td>C1</td><td>11.4 ± 0.7</td><td>8.71 <sup>+0.10</sup><sub>-0.08</sub></td><td>0.2 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.6 <sup>+0.1</sup><sub>-0.2</sub></td><td>0.38 <sup>+0.21</sup><sub>-0.20</sub></td></tr><tr><td>C2</td><td>9.2 ± 0.9</td><td>8.02 <sup>+0.15</sup><sub>-0.43</sub></td><td>0.3 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.5 <sup>+0.1</sup><sub>-0.4</sub></td><td>2.86 <sup>+1.37</sup><sub>-2.99</sub></td></tr><tr><td>C3</td><td>9.8 ± 0.8</td><td>8.55 <sup>+0.09</sup><sub>-0.15</sub></td><td>0.2 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.5 <sup>+0.1</sup><sub>-0.3</sub></td><td>0.56 <sup>+0.30</sup><sub>-0.34</sub></td></tr><tr><td>Smooth</td><td>58.7 ± 7.3</td><td>9.35 <sup>+0.07</sup><sub>-0.05</sub></td><td>0.4 <sup>+0.4</sup><sub>-0.2</sub></td><td>0.7 <sup>+0.1</sup><sub>-0.2</sub></td><td>0.17 <sup>+0.18</sup><sub>-0.09</sub></td></tr><tr><td>Integrated</td><td>89.2 ± 8.3</td><td>9.36 <sup>+0.09</sup><sub>-0.11</sub></td><td>1.3 <sup>+1.0</sup><sub>-0.4</sub></td><td>0.4 <sup>+0.2</sup><sub>-0.1</sub></td><td>0.56 <sup>+0.45</sup><sub>-0.22</sub></td></tr><tr><td>JADES-JEMS-6885</td><td>3.698</td><td>C1</td><td>39.9 ± 1.6</td><td>9.19 <sup>+0.07</sup><sub>-0.03</sub></td><td>1.5 <sup>+0.1</sup><sub>-0.2</sub></td><td>0.6 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.96 <sup>+0.16</sup><sub>-0.14</sub></td></tr><tr><td></td><td></td><td>C2</td><td>16.8 ± 0.8</td><td>8.30 <sup>+0.10</sup><sub>-0.15</sub></td><td>1.2 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.3 <sup>+0.1</sup><sub>-0.1</sub></td><td>6.01 <sup>+1.47</sup><sub>-2.13</sub></td></tr><tr><td></td><td></td><td>C3</td><td>12.2 ± 0.8</td><td>8.13 <sup>+0.08</sup><sub>-0.09</sub></td><td>1.2 <sup>+0.1</sup><sub>-0.2</sub></td><td>0.1 <sup>+0.1</sup><sub>-0.1</sub></td><td>8.89 <sup>+1.79</sup><sub>-2.36</sub></td></tr><tr><td></td><td></td><td>Smooth</td><td>169.0 ± 13.0</td><td>9.62 <sup>+0.07</sup><sub>-0.21</sub></td><td>3.6 <sup>+1.1</sup><sub>-0.6</sub></td><td>0.4 <sup>+0.2</sup><sub>-0.2</sub></td><td>0.86 <sup>+0.29</sup><sub>-0.14</sub></td></tr><tr><td></td><td></td><td>Integrated</td><td>238.1 ± 14.0</td><td>9.87 <sup>+0.03</sup><sub>-0.07</sub></td><td>7.4 <sup>+1.1</sup><sub>-0.9</sub></td><td>0.4 <sup>+0.2</sup><sub>-0.2</sub></td><td>0.99 <sup>+0.16</sup><sub>-0.20</sub></td></tr><tr><td rowspan="4">JADES-JEMS-16296</td><td rowspan="4">3.53</td><td>C1</td><td>42.5 ± 1.3</td><td>8.82 <sup>+0.02</sup><sub>-0.03</sub></td><td>2.7 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.5 <sup>+0.1</sup><sub>-0.1</sub></td><td>4.08 <sup>+0.24</sup><sub>-0.32</sub></td></tr><tr><td>C2</td><td>7.1 ± 0.9</td><td>7.88 <sup>+0.05</sup><sub>-0.08</sub></td><td>0.4 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.1 <sup>+0.1</sup><sub>-0.1</sub></td><td>5.27 <sup>+1.45</sup><sub>-1.63</sub></td></tr><tr><td>Smooth</td><td>80.8 ± 7.3</td><td>9.37 <sup>+0.02</sup><sub>-0.01</sub></td><td>0.9 <sup>+0.1</sup><sub>-0.2</sub></td><td>0.9 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.38 <sup>+0.04</sup><sub>-0.08</sub></td></tr><tr><td>Integrated</td><td>130.6 ± 8.4</td><td>9.65 <sup>+0.05</sup><sub>-0.06</sub></td><td>5.4 <sup>+1.2</sup><sub>-0.9</sub></td><td>0.5 <sup>+0.3</sup><sub>-0.2</sub></td><td>1.20 <sup>+0.30</sup><sub>-0.26</sub></td></tr><tr><td rowspan="4">JADES-JEMS-14436</td><td rowspan="4">3.67</td><td>C1</td><td>87.8 ± 1.8</td><td>9.58 <sup>+0.05</sup><sub>-0.02</sub></td><td>3.9 <sup>+0.2</sup><sub>-0.2</sub></td><td>0.4 <sup>+0.1</sup><sub>-0.1</sub></td><td>1.02 <sup>+0.12</sup><sub>-0.07</sub></td></tr><tr><td>C2</td><td>38.5 ± 0.9</td><td>9.32 <sup>+0.11</sup><sub>-0.06</sub></td><td>5.8 <sup>+0.6</sup><sub>-0.7</sub></td><td>0.3 <sup>+0.2</sup><sub>-0.1</sub></td><td>2.77 <sup>+0.75</sup><sub>-0.50</sub></td></tr><tr><td>Smooth</td><td>169.9 ± 17.0</td><td>10.05 <sup>+0.04</sup><sub>-0.06</sub></td><td>2.8 <sup>+0.6</sup><sub>-0.6</sub></td><td>0.9 <sup>+0.1</sup><sub>-0.2</sub></td><td>0.24 <sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td>Integrated</td><td>296.3 ± 18.0</td><td>10.34 <sup>+0.07</sup><sub>-0.09</sub></td><td>11.8 <sup>+2.5</sup><sub>-1.9</sub></td><td>0.8 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.53 <sup>+0.14</sup><sub>-0.14</sub></td></tr><tr><td rowspan="4">JADES-JEMS-11059</td><td rowspan="4">4.07</td><td>C1</td><td>47.8 ± 1.2</td><td>9.32 <sup>+0.02</sup><sub>-0.03</sub></td><td>2.0 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.7 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.95 <sup>+0.06</sup><sub>-0.08</sub></td></tr><tr><td>C2</td><td>24.0 ± 1.1</td><td>7.97 <sup>+0.05</sup><sub>-0.07</sub></td><td>1.2 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.4 <sup>+0.1</sup><sub>-0.1</sub></td><td>12.85 <sup>+1.82</sup><sub>-2.33</sub></td></tr><tr><td>C3</td><td>8.2 ± 0.5</td><td>8.22 <sup>+0.22</sup><sub>-0.39</sub></td><td>0.4 <sup>+0.1</sup><sub>-0.1</sub></td><td>0.5 <sup>+0.1</sup><sub>-0.3</sub></td><td>2.41 <sup>+1.36</sup><sub>-2.24</sub></td></tr><tr><td>Smooth</td><td>164.0 ± 11.6</td><td>9.81 <sup>+0.07</sup><sub>-0.12</sub></td><td>7.3 <sup>+1.7</sup><sub>-2.5</sub></td><td>0.3 <sup>+0.4</sup><sub>-0.2</sub></td><td>1.13 <sup>+0.32</sup><sub>-0.49</sub></td></tr><tr><td></td><td></td><td>Integrated</td><td>244.2 ± 12.9</td><td>9.59 <sup>+0.08</sup><sub>-0.07</sub></td><td>9.6 <sup>+0.7</sup><sub>-0.6</sub></td><td>0.4 <sup>+0.1</sup><sub>-0.2</sub></td><td>2.46 <sup>+0.48</sup><sub>-0.42</sub></td></tr></tbody></table>

**Notes.**

<sup>a</sup> Fluxes from aperture photometry in PSF-matched NIRCam/F150W images.

<sup>b</sup> Reported values are from the PROSPECTOR SED fitting with our fiducial model (Section 2.4).

<sup>c</sup> Instantaneous SFR from SED fitting, i.e., SFR in the first lookback time bin (30 Myr) of a nonparametric SFH.

the smooth component surrounding the clumpy substructures in each one of the galaxies. Ideally, in the FORCEPHO fitting, we could add an additional morphological component to model this smooth component together with the clumpy substructures. However, we find that a single Sérsic profile is unable to describe the smooth component, because its morphology is rather complex with large-scale asymmetries. In principle, we could assume a more complex model for it, e.g., a summation of several Sérsic profiles. But such a model is purely arbitrary, and it would be hard to understand and control the systematic errors behind it. We therefore decide not to model the smooth component in the FORCEPHO fitting. One concern is then how this can affect the morphological fits to the clumpy substructures. While it can be significant for faint clumps C2 and C3, we argue that the effect should be minor for C1 (which is the focus of this work) because the ratio of the median residual flux (primarily from the smooth component) to the peak flux of C1 is only ~0.1%–4%.

**3. Results**

In this Section, we present the spatially resolved stellar populations of the six galaxies. In particular, we focus on the properties of the red, massive core C1 in each one of the galaxies, and discuss the constraints on its formation mechanisms.

**3.1. Spatially Resolved Stellar Populations**

Table 2 presents the physical properties derived from the fiducial PROSPECTOR model described in Section 2.4. To check the SED fitting results, we first compare the stellar mass and SFR derived using the integrated flux with the sum of that measured for individual components, i.e., substructures plus the smooth component. We note that these two values are not necessarily consistent with each other, because summation of the measurements of *N* individual components is equivalent to fitting a galaxy’s SED with (N − 1)× more free parameters
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Spatially resolved stellar populations of the three galaxies with $z_{\text{spec}}$. Left: the RGB images are the same as those shown in Figure 3, except that here we zoom in to the central $1''/3 \times 1''/3$ region. The middle and right panels show the SED fitting results from our fiducial PROSPECTOR model (Section 2.4). Middle: the normalized (i.e., $F_{\nu}^{\text{F150W}} = 1$) photometry and best-fit SEDs of individual structural components defined in Section 2.3. These include best-fit SEDs for (1) visually identified clumpy substructures, which are color-coded accordingly based on the colors used in the RGB images for individual apertures, (2) the smooth component (orange), and (3) integrated flux (black). The inset highlights the wavelength range around $2\mu\text{m}$ covered by NIRCam F150W, F182M, F200W, and F210M, and HST/WFC3 F160W filters. The magenta shaded region marks the rest-frame wavelength range between 3645 and 4100 Å, which contains key spectral features sensitive to stellar age, such as the Balmer jump and D4000 break (A.G. Bruzual 1983; M. L. Balogh et al. 1999). Right: reconstructed nonparametric SFHs of individual substructures. Each SFH is normalized with its peak SFR, i.e., normalized SFR(t) = SFR(t)/max(SFR(t)), to better show the difference in the shape of SFHs. In each subpanel, we overplot the best-fit nonparametric SFHs for the smooth component (orange) and integrated flux (black). For better illustration, here we only plot the best-fit SFHs with solid lines, while in Appendix C, we also plot the $1\sigma$ uncertainties.</figcaption></figure>

than using its integrated flux alone. Yet, excellent agreement between the two measurements is seen in Figure 5, showing the consistency and robustness of our SED fitting. In Appendix B, we further compare the coadded SFH of individual structural components, with the SFH obtained from SED fitting of integrated flux. The two SFHs are consistent with each other within uncertainty.

Figures 6 and 7 show the spatially resolved SED fitting results in detail. The SED shapes of individual structural components differ across the entire rest-frame UV to NIR wavelengths, suggesting that stellar populations vary significantly within each one of the galaxies. This highlights the importance of having spatially resolved measures of stellar populations to get a comprehensive view on the mass assembly processes in high-redshift galaxies.

We refer the readers to Appendix D for detailed descriptions of SED fitting results of individual galaxies. In the list below, we only summarize the key findings. In each one of the galaxies,

1. the red stellar core C1 is the dominant, most massive one among all identified clumpy substructures. Relative to the off-center, minor clumps C2 and C3, C1 contains older stellar populations, spanning a range in mass-weighted stellar age of 0.4–0.7 Gyr, and has a lower specific SFR. In all galaxies, C1 has a common feature in its reconstructed SFH—it experienced a recent major star formation episode that started ~0.5–1 Gyr and peaked at ~0.1–0.3 Gyr prior to the time of observation.

2. compared to C1, the minor, off-center clumps C2 and C3 are less massive by $\gtrsim 0.5$ dex. They have bluer colors and larger specific SFRs, and feature a generally rising SFH.
<figure><img src="image_8.png" /><figcaption>Figure 7. Similar to Figure 6, but for the three galaxies with $z_{phot}$ shown in Figure 4. Note that the C2 of JADES-JEMS-11059 is not covered by its segmentation map (the bottom row of Figure 4), but we still include it in the analysis because its photometric redshift is consistent with that of the central galaxy.</figcaption></figure>

3. the stellar population of the smooth component has a comparable or older stellar age than C1. The smooth component contributes the most stellar mass and SFR to the entire galaxy, but it has smaller surface densities (i.e., mass/SFR per area) compared to the clumpy substructures.

## 3.2. Stellar Masses of the Substructures

The stellar mass of clumpy substructures contains key information on their physical origins (e.g., A. Dekel et al. 2009; R. Genzel et al. 2011; Y. Guo et al. 2012; F. Bournaud et al. 2014; Y. Guo et al. 2015; M. Dessauges-Zavadsky & A. Adamo 2018). Directly comparing their stellar masses with the mass threshold of gravitational collapse due to local hydrodynamical instabilities ($M_{Toomre}$; A. Toomre 1964, see Section 3.3 for details) will immediately tell us that if their formation is consistent with in situ fragmentation in gaseous disks.

Accurately measuring the stellar mass of clumps, however, is very challenging, because they are embedded in a smooth component, e.g., a galaxy disk. We thus need to properly subtract the contribution from the smooth component at the

clump locations; otherwise, the stellar-mass measures will be biased (M. Huertas-Company et al. 2020). In our case, this can be formularized as

$$M_{C}^{\text{int}} = \eta \cdot (M_{C} - M_{S}^{C}), \qquad (4)$$

where $M_{C}^{\text{int}}$ is the intrinsic (i.e., true) stellar mass of a clump; $M_{C}$ is the stellar mass of that clump derived from SED fitting with fluxes directly from aperture photometry (Section 2.3 and Section 2.4); $M_{S}^{C}$ is the contribution from the smooth component to the stellar-mass measure at the clump location; and finally $\eta$ accounts for aperture correction.

Reliably estimating $M_{S}^{C}$ and $\eta$ requires information on intrinsic shapes of clumps and the smooth component, which unfortunately are not known a priori, and it becomes increasingly challenging when clumps are fainter. However, we reiterate that the primary goal of this work is to understand the formation of massive stellar cores (i.e., C1), rather than to, e.g., measure the stellar-mass function of clumps. If the stellar mass of a clump is less than the threshold mass even before subtracting $M_{S}^{C}$, i.e., $\eta M_{C} < M_{\text{Toomre}}$, its formation will then be fully consistent with in situ fragmentation. Robustly determining the stellar masses for these low-mass, faint clumps (mostly C2 and C3) is beyond the scope of this paper, and in such cases
we simply use $\eta M_C$ as a strict upper limit (Figure 9 and Section 3.3).

In the remainder of this Section, the main focus will be on estimating intrinsic stellar masses of clumps (C1 in particular) more massive than $M_{\text{Toomre}}$. We estimate them using two different empirical methods (Sections 3.2.1 and 3.2.2), and with different aperture correction factors (Section 3.2.3). In short, our conclusion, that intrinsic stellar masses of C1 are larger than $M_{\text{Toomre}}$, does not depend on the choice of method.

### 3.2.1. Estimating $M_C^{\text{int}}$ by Directly Subtracting the Mass Contribution from the Smooth Component

As we showed earlier in Figure 5, because the coadded measurements of stellar mass and SFR are consistent with those derived using integrated flux, this supports our decision to perform the direct subtraction in mass. If we assume that the mass-to-light ($M/L$) ratio is constant and stellar mass is uniformly distributed within the smooth component, then $M_S^C$ simply is the stellar mass of the smooth component derived via SED fitting using fluxes directly from aperture photometry ($M_S$) and scaled by the area ratio of the clump ($A_C$) to the smooth component ($A_S$). The intrinsic stellar mass of the clump is then

$$M_C^{\text{int}} = \eta \cdot \left( M_C - \frac{A_C}{A_S} M_S \right). \qquad (5)$$

Because stellar masses of individual components were derived by SED fitting (Section 2.4), the biggest advantage of this approach is that the spatial variation in $M/L$ has already been largely taken into account, despite not being on a pixel-by-pixel basis. In Figure 9, we plot the resulting $M_C^{\text{int}}$ as filled symbols with error bars. For all C1, we find $M_C^{\text{int}}/M_{\text{Toomre}} \gtrsim 2$, arguing against formation via in situ fragmentation. We use the intrinsic stellar masses of C1 derived from this approach as the fiducial ones for the discussion in Section 4.

### 3.2.2. Estimating $M_C^{\text{int}}$ by Subtracting the Flux Contribution from the Smooth Component

Differing from directly subtracting mass, alternatively we can first subtract the flux contribution from the smooth component to the clump location, and then estimate $M_C^{\text{int}}$. The Equation (4) thus becomes

$$M_C^{\text{int}} = \eta \cdot \left( M_C - \frac{A_C \cdot f_S^C}{F_C} M_C \right), \qquad (6)$$

where $f_S^C$ is the flux contribution per pixel from the smooth component to the clump location, and $F_C$ is the flux of the clump from aperture photometry.

For each one of the galaxies, we estimate $f_S^C$ by subtracting from its F210M image the best-fit FORCEPHO models (Section 2.5) of all clumpy substructures. We then do a $3\sigma$ clipping on the residual image and use the median residual pixel value as $f_S^C$. Because we ignore the smooth component during the FORCEPHO fitting (Section 2.5), instead of using the residual image, we also estimate $f_S^C$ in a more aggressive way. We first mask out pixels within $r=2 \times R_e$ circular regions of individual clumpy substructures from the galaxy's segmentation map, and then estimate $f_S^C$ as the median pixel value of the masked image (i.e., without subtracting the best-fit

FORCEPHO models of clumpy substructures). This latter estimate certainly overestimates $f_S^C$ because the flux contribution from clumps to the smooth component is not fully removed. Even so, we find that $f_S^C$ only increases by ~3%–7%, which is small and will not cause any substantial changes in our results.

In Figure 9, we use hollow symbols with dashed edges to show the $M_C^{\text{int}}$ of C1 from this approach. The $M_C^{\text{int}}$ from this approach are larger than the fiducial ones, which is expected because Equation (6) assumes the same $M/L$ for C1 and the smooth component. As Figure 6 and 7 show, in each one of the galaxies, C1 has a redder rest-frame (B-V) color, as probed by F182M-F277W, than the smooth component. Because $M/L$ increases as (B-V) becomes redder (e.g., S. S. McGaugh & J. M. Schombert 2014), using Equation (6) can lead to an under-subtraction of $M_S^C$. This issue can be mitigated by doing the FORCEPHO analysis across all filters, as opposed to only using the F210M filter, and then taking the $M/L$ variation into account with the FORCEPHO-derived model photometry. Following this, we find that the intrinsic masses of C1 can change by <20%, a sufficiently small amount that does not substantially change any of our results.<fnref n="26" />

### 3.2.3. Estimating the Aperture Correction Factor $\mu$

By default, we use the best-fit Sérsic profile from the FORCEPHO fitting (Section 2.5) to derive the aperture correction factor $\mu$. As the posterior plots of Figure 8 show, the Sérsic profiles of C1 are tightly constrained, with typical uncertainties of Sérsic index $n$ and effective radius $R_e$ of <10%. C1 are resolved/marginally resolved in NIRCam/F210M images (PSF FWHM ≈0.07), spanning a range in $R_e$ of 0.07–0.14, corresponding to 0.5–1 kpc at the median redshift $z=3.7$. For minor clumps C2 and C3, their Sérsic profiles are much less constrained. We also check our results by adopting a rather conservative aperture correction, namely, assuming the clumps are point sources, which we already knew does not accurately describe C1's morphologies. Using the point-source aperture correction can lead to a decrease in stellar mass by a factor of ~1.5–2.1, which, still, cannot substantially change our conclusions because all C1 have $M_C^{\text{int}}/M_{\text{Toomre}} \gtrsim 2$ (see Figure 9 and Section 3.3 below).

## 3.3. Origins of the Substructures

A prominent feature of high-redshift galaxies is the presence of giant stellar clumps that are far more massive than individual star clusters observed in nearby galaxies (e.g., D. M. Elmegreen et al. 2007; R. Genzel et al. 2008, 2011; S. Wuyts et al. 2012; Y. Guo et al. 2012, 2018). One scenario to form such giant clumpy substructures is through in situ fragmentation due to VDIs (e.g., M. Noguchi 1999; A. Immeli et al. 2004; F. Bournaud et al. 2007; B. G. Elmegreen et al. 2008; A. Dekel et al. 2009; D. Ceverino et al. 2010; S. Inoue et al. 2016). In a rotating gaseous disk, gravity needs to overcome both gas pressure and

<footnote n="26">²⁶ We could use the FORCEPHO model photometry from the outset to do SED fitting, instead of using aperture photometry. But we argue that using the model photometry would then make the accuracy of the measured SED sensitive to the assumption of the intrinsic shapes of the clumps, which could deviate from a single Sérsic profile and may also depend on wavelength. These systematics are difficult to control, and can be particularly significant for the inferences of stellar age and SFH. Therefore, in this work we decide to stay empirical and use the properties of stellar populations measured with fluxes from aperture photometry.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Morphological analysis for the individual substructures. Each substructure is modeled as a Sérsic profile using the code FORCEPHO that adopts an MCMC framework to perform PSF-convolved analysis in individual exposures to better estimate the parameter uncertainties (Section 3.2). The best-fit centroids of individual substructures are marked with the star symbols. The left three columns show the F210M cutout, best-fit model, and residual, respectively. The fourth column shows the residual χ map, which is defined as $χ = (data − model)/error$, from which we clearly see the smooth component surrounding the clumpy substructures. The last column shows the posterior corner plot of Sérsic index *n* vs. circularized effective (i.e., half-light) radius *R*<sub>e</sub>.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The scatter plot of the intrinsic stellar masses of clumpy substructures (y-axis) vs. the total stellar masses of their host galaxies (x-axis). The black solid line marks the prediction of Toomre mass $M_{\text{Toomre}}$, which represents the largest stellar mass that a clump can form through Toomre instability in an unstable, gaseous disk (see Section 3.3). The dashed and dotted lines mark $2 \times$ and $3 \times M_{\text{Toomre}}$, respectively. The upper-left gray shaded region marks the unphysical part where $M_{\text{C}}^{\text{int}} > M_{*}$. The symbols used for individual galaxies, as well as the corresponding color scheme used for individual clumpy substructures, are labeled on the right. For low-mass clumps whose stellar masses are already less than $M_{\text{Toomre}}$ before subtracting the contribution from the surrounding smooth component (Section 3.2), we plot the upper limits of their $M_{\text{C}}^{\text{int}}$. Similar to previous plots, the stellar cores, C1, are color-coded in red. The filled symbols with error bars show the $M_{\text{C}}^{\text{int}}$ estimated via the method of Section 3.2.1, while the hollow symbols with dashed edges show the $M_{\text{C}}^{\text{int}}$ estimated via a different method described in of Section 3.2.2. We find that the intrinsic stellar masses of C1 are significantly larger than $M_{\text{Toomre}}$ in all six galaxies, suggesting that C1 may not form through in situ fragmentation due to Toomre instability.</figcaption></figure>

shear forces produced by the relative motion between fluid parcels. A. Toomre (1964) showed that in situ fragmentation can happen if the parameter $Q$ is less than unity (Section 1). Consequently, clumps have a characteristic mass that represents the largest stellar mass that a clump can form via Toomre instability, a.k.a. the Toomre mass. Following R. Genzel et al. (2011), the Toomre mass can be calculated as

$$M_{\text{Toomre}} \approx 5 \times 10^9 \left( \frac{f_{\text{young}}}{0.4} \right)^2 \left( \frac{M_{\text{disk}}}{10^{11} M_{\odot}} \right) M_{\odot}, \quad (7)$$

where $f_{\text{young}}$ is the mass fraction of components forming stars, and $M_{\text{disk}}$ is the total baryonic mass of a disk.

To estimate $M_{\text{Toomre}}$ for our sample, we assume $f_{\text{young}} \approx f_{\text{gas}}$, where $f_{\text{gas}}$ is the gas-to-baryonic mass fraction, i.e., $M_{\text{gas}}/(M_{\text{gas}} + M_{*})$. Strictly speaking, $M_{\text{gas}}$ here should be the summation of molecular and atomic gas masses in galaxy disks, but we assume $M_{\text{gas}} \approx M_{\text{molgas}}$ for the following reason. While directly measuring atomic gas in galaxies beyond the local Universe is technically challenging, studies using the damped Ly$\alpha$ absorbers toward high-redshift quasars infer a generally very weak, if any, evolution of the cosmic atomic hydrogen density ($\rho_{\text{H}}$) up to $z \sim 5$ (e.g., S. M. Rao et al. 2006; J. X. Prochaska & A. M. Wolfe 2009; C. Péroux & J. C. Howk 2020; F. Walter et al. 2020). Meanwhile, observations with different probes of molecular gas, including CO, far-IR SED, and submillimeter dust continuum, have reached a consensus, at least up to $z = 4$, that the molecular-gas mass of galaxies rapidly evolves at a rate of $(1+z)^{3 \pm 1}$ (e.g., R. Genzel et al. 2015; N. Scoville et al. 2017; L. J. Tacconi et al. 2018; D. Liu et al. 2019). Putting together the weak $\rho_{\text{H}}$ and strong molecular-gas evolutions, because the mass ratio of atomic gas divided by molecular gas is $\approx 2-3$ at $z \sim 0$ (e.g., A. Saintonge et al. 2011; B. Catinella et al. 2013), it is reasonable to assume $M_{\text{gas}} \approx M_{\text{molgas}}$ (hence $f_{\text{gas}} \approx f_{\text{molgas}}$) in

galaxies at $z \gtrsim 1$ (see a similar argument from L. J. Tacconi et al. 2018 for this assumption).

We estimate $f_{\text{molgas}}$ for the six galaxies using the best-fit relationship from L. J. Tacconi et al. (2020; see their Table 3), who combined a number of surveys of galaxies’ molecular gas content across cosmic time, and also included the dependence of the relationship on stellar mass and the distance to the star-forming main sequence. We note, however, that almost all of those high-redshift surveys were designed for galaxies with $M_{*} > 10^{10} M_{\odot}$, either on or above the star-forming main sequence. Because five out of the six galaxies presented in this work have $M_{*} \approx 10^{9.3-10} M_{\odot}$ (Table 1), we use the $f_{\text{molgas}}$ derived for $M_{*} = 10^{10} M_{\odot}$ galaxies. As we will show in Section 4.2, all six galaxies are below the star-forming main sequence. Using their distances to the star-forming main sequence, we find $f_{\text{gas}} \approx 0.5$ for our sample, similar to those observed in star-forming galaxies on the main sequence at $z \sim 2$ (e.g., L. J. Tacconi et al. 2008; E. Daddi et al. 2010; L. J. Tacconi et al. 2010). We note that $f_{\text{molgas}}$ increases with decreasing stellar mass. If the scaling relationship of $f_{\text{molgas}}$ from L. J. Tacconi et al. (2020) can be safely used for galaxies with $M_{*} < 10^{10} M_{\odot}$, we would get $f_{\text{molgas}} \sim 0.6$ assuming $M_{*} = 10^{9.3} M_{\odot}$, the lowest stellar mass of the six galaxies. This would increase $M_{\text{Toomre}}$ by a factor of $(0.6/0.5)^2 = 1.44$, which however is small compared to the stellar mass of C1 and hence does not change our conclusion about the formation of C1 detailed below.

With $f_{\text{gas}}$ in hand, we then get $M_{\text{disk}} = M_{*} + M_{\text{gas}} = M_{*}/f_{\text{gas}}$, and finally use Equation (7) to calculate $M_{\text{Toomre}}$. In Figure 9, we plot the relationship between $M_{\text{Toomre}}$ and $M_{*}$ as the black solid line. In what follows, we discuss possible physical origins of the clumpy substructures.

To begin, Figure 9 shows that the most minor clumps (C2 and C3) are less massive than $M_{\text{Toomre}$; hence, they are consistent with in situ formation through gravitational
fragmentation in gaseous, turbulent disks. These clumps are bluer, and have enhanced specific star formation rates (sSFRs) compared to C1 and the smooth component (Table 2). These properties are similar to the star-forming clumps found in galaxies at lower redshifts $1 < z < 3$ from earlier HST studies (e.g., B. G. Elmegreen & D. M. Elmegreen 2005; D. M. Elmegreen et al. 2007; S. Wuyts et al. 2012). Quantitatively, these individual minor clumps contribute 5%–50% of the total SFR, and <10% of the total stellar mass of their host galaxies, in broad agreement with earlier studies of UV-bright clumps at $1 < z < 3$ (e.g., Y. Guo et al. 2012). In addition, they have overall rising SFHs that started ≈0.3 Gyr ago, which is also consistent with the young stellar ages found in $2 < z < 3$ star-forming clumps (D. M. Elmegreen et al. 2009; N. M. Förster Schreiber et al. 2011; R. Genzel et al. 2011; Y. Guo et al. 2018).

Apart from C1, two other clumpy substructures are more massive than $M_{\text{Toomre}}$, i.e., the C3 of JADES-JEMS-15157 (blue X in Figure 9) and the C2 of JADES-JEMS-14436 (green hexagon in Figure 9). Regarding the former, its intrinsic stellar mass is $\sim 1.5 \times M_{\text{Toomre}}$. We notice, however, that it might contain two distinct clumps, which is particularly clear in the F210M image but becomes less obvious in other filters (Figure 3). After taking this into account, its formation is then still consistent with in situ fragmentation. Regarding the latter, i.e., the C2 of JADES-JEMS-14436, we see a clear tidal feature associated with this substructure in the NIRCam SW's filter images (Figure 4 and Appendix D). Therefore, the physical origin of this substructure is very likely related to a minor merger with a mass ratio of 1:10 (Table 2). The minor mergers can have two broad effects on the clump formation. First, minor mergers can disturb the gas distribution in galaxy disks, trigger local starbursts, and hence induce clump formation (F. Bournaud et al. 2014; S. Inoue et al. 2016). Second, the clumps can be minor mergers themselves, i.e., the clumps have an ex situ origin. Using cosmological hydrodynamical simulations, N. Mandelker et al. (2014) found that the ex situ clumps have a typical stellar age of 0.5–3 Gyr and sSFR of $0.1–2 \text{ Gyr}^{-1}$, which are in quantitative agreement with the C2 of JADES-JEMS-14436 (Figure 7 and Appendix D).

Finally, at the core of this study is the physical origin of C1—the stellar core. In each one of the six galaxies, C1 is redder than the other parts of the galaxy, and has significantly larger stellar mass than $M_{\text{Toomre}}$, with a typical ratio of $M_C^{\text{int}}/M_{\text{Toomre}} \gtrsim 2$. This excessive stellar mass is strong evidence that C1 may not formed via a single in situ fragmentation in gaseous disks. Relative to those minor clumps discussed above, C1 has a lower sSFR and is closer to the light centroid of NIRCam images at $\sim 4 \mu\text{m}$. All of these are consistent with C1 being a forming stellar core/(pseudo-)bulge, which will be discussed in detail in the remainder of this paper.

# 4. Discussion

## 4.1. The Formation of C1: Early (Pseudo-)Bulge Formation?

In the local Universe, bulges and pseudo-bulges are believed to have different origins, with the former primarily formed through major mergers (e.g., L. Hernquist 1989; J. E. Barnes & L. Hernquist 1996; P. F. Hopkins et al. 2009) and the latter formed through secular evolution in galaxy disks (e.g., V. P. Debattista et al. 2004; T. Okamoto 2013; E. Athanassoula et al. 2015). According to the classification for $z \sim 0$ galaxies,

the light distribution of a classic bulge is similar to a spheroid ($n=4$), while it becomes flattened (smaller $n$) for a pseudo-bulge (J. Kormendy & R.C.J. Kennicutt 2004). Therefore, from a purely morphological perspective (Figure 8), the C1 of JADES-JEMS-13396 ($n>2$) is more like a classic bulge, while in the other five galaxies, C1 are consistent with pseudo-bulges; although, we stress that without dynamical measures it remains difficult to definitively distinguish the two.

At high redshifts, early numerical simulations showed that stellar cores formed through the coalescence of inspiral stellar clumps should have $n \sim 2–4$, which is consistent with classic bulges (B. G. Elmegreen et al. 2008). Later simulations however, after taking into account the gas kinematics of individual clumps, showed that the remnant of the clump coalescence should still have significant rotation, arguing that the high-redshift cores are more consistent with being pseudo-bulges (D. Ceverino et al. 2012; S. Inoue & T. R. Saitoh 2012). Regardless, those simulations consistently showed that the coalescence of clumps is an efficient way to grow dense central cores in gas-rich galaxies at high redshifts. Indeed, the coexistence of C1 with other minor, off-center star-forming clumps (C2 and C3) in the six galaxies is consistent with this picture.

The reconstructed SFHs provide another key constraint on the physical origin of C1. If they form through the coalescence of inspiral clumps, the timescale for the clumps' inward migration will be ≈0.5 Gyr (e.g., M. Noguchi 1999; A. Immeli et al. 2004; F. Bournaud et al. 2007; R. Genzel et al. 2008, 2011). This is consistent with the SED fitting that the mass-weighted stellar ages of C1 span a range of 0.4–0.7 Gyr (Section 3.1). The clump migration scenario also predicts the radial age gradient of clumps. Indeed, those off-center minor clumps (C2 and C3, if present) have significantly younger SFHs (Figures 6 and 7), spanning a range in mass-weighted stellar age of 0.1–0.5 Gyr (Table 2). These are in quantitative agreement with the scenario of growing the stellar core of massive galaxies at high redshift through VDI, i.e., one of the processes associated with gas-rich compaction initially proposed by A. Dekel et al. (2009); also see Section 1 for a detailed description), which predicts that clumps should have stellar ages <0.5 Gyr at the outskirts and be older at the center.

One important question about the efficacy of the inward migration scenario is the uncertain lifetime of clumps. Can clumps formed in disks survive stellar feedback until they coalesce at the center? Some simulations with intense stellar feedback models showed that clumps are short-lived, and should be self-destroyed within <50–100 Myr (P. F. Hopkins et al. 2012; S. Genel et al. 2012). In contrast, other simulations with different stellar feedback models showed that they can produce relatively longer-lived stellar clumps with a typical age of a few hundred megayears (J. Perez et al. 2013; D. Ceverino et al. 2014; F. Bournaud et al. 2014). The stellar ages of 0.1–0.5 Gyr measured for the minor clumps C2 and C3 favor a relatively longer lifetime of clumps than the predictions from intense stellar feedback models. Similar conclusions were also reached by earlier studies of $z \sim 2$ star-forming clumps (S. Wuyts et al. 2012; Y. Guo et al. 2012). Arguably, however, the measured stellar ages of C2 and C3 still seem to be shorter than the ≈0.5 Gyr required for the inward migration. As pointed out by F. Bournaud (2016), however, the observed stellar ages of clumps are only lower limits to their true ages, because during their evolution, clumps continuously exchange
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The star-forming main sequence of $3 < z < 4.5$ galaxies from JEMS. The stellar mass and SFR are derived from SED fitting with integrated flux. Each galaxy is color-coded using its best-fit photometric redshift. The green solid line shows the star-forming main sequence at $z = 3.7$ for $M_{*} > 10^{8.5} M_{\odot}$ galaxies from P. Popesso et al. (2023). We extrapolate it to $M_{*} < 10^{8.5} M_{\odot}$, and plot it as the green dashed line. The green shaded region marks the ±0.3 dex range. Also plotted as the magenta solid line is the star-forming main sequence at $z = 2.7$ for $M_{*} > 10^{10.2} M_{\odot}$ galaxies from J. Leja et al. (2022), who also used PROSPECTOR with a similar SED model of this work. Similarly, we also extrapolate it to $M_{*} < 10^{10.2} M_{\odot}$, and plot it as the magenta dashed line. The six galaxies presented in this work are marked with golden hexagons. They are below the star-forming main sequence by ≈0.2–0.7 dex.</figcaption></figure>

mass with the surrounding ISM through outflows and inflows (F. Bournaud et al. 2007; A. Dekel & M. R. Krumholz 2013; F. Bournaud et al. 2014; V. Perret et al. 2014). According to the zoom-in simulations from F. Bournaud et al. (2014), a clump with an observed stellar age of 0.1–0.2 Gyr corresponds to a true age of 0.3–0.5 Gyr. Taking this into account, the inferred stellar populations for the minor clumps thus are still in line with the scenario where the central regions of high-redshift galaxies are built through the inspiral of clumps.

We also note that gas-rich major mergers (another process associated with gas-rich compaction), if present, can help efficiently grow central cores/bulges at high redshifts as well. However, we do not see clear evidence of major mergers in the six galaxies presented here. The effect induced by recent minor mergers, e.g., the formation of bars, can also contribute to the buildup of (pseudo-)bulges (e.g., F. Bournaud et al. 2005; M. C. Eliche-Moral et al. 2011; J. Guedes et al. 2013). Given the complex morphologies of the galaxies, however, it is hard to distinguish an in situ stellar clump from a merging low-mass galaxy. Minor mergers thus might also contribute to the development of central regions in the six galaxies.

Finally, we notice that C1 has a comparable/younger mass-weighted stellar age than the smooth component in each one of the galaxies (Table 2). As the reconstructed SFHs show (see Figures 6 and 7), this is because the smooth component started its major mass assembly either at about the same time as C1, or earlier, i.e., >1 Gyr prior to the time of observation. Interestingly, recent studies of galactic archeology found that the Milky Way’s bulge is predominately composed of $\gtrsim 9-10$ Gyr-old stars (corresponding to a formation redshift $z \lesssim 2-3$; e.g., B. Barbuy et al. 2018; S. Hasselquist et al. 2020), while the formation of thick-disk stars likely started as early as at $z > 6$, i.e., 13 Gyr ago (M. Xiang & H.-W. Rix 2022). These are seemingly similar to the reconstructed SFHs of C1 and the smooth component. It is thus possible that we are witnessing the early formation of structural components that will

eventually evolve into the bulge and thick-disk components in the Milky Way–like galaxies at $z \sim 0$; although we caution that the six galaxies are likely hosted by dark matter halos ≈0.5–1.5 dex more massive<fnref n="27" /> than the Milky Way’s.

## 4.2. Implications for Galaxy Quenching

We now advance the discussion by connecting what we found about the six galaxies to the general picture of galaxy quenching.

To begin, in Figure 10 we present the star-forming main sequence, i.e., SFR versus stellar mass, for all 301 galaxies at $3 < z < 4.5$ in our parent sample from JEMS (Section 2.1). The SFR and stellar-mass measures shown in the plot are from our PROSPECTOR SED fitting with the same model as that assumed for the six galaxies (Section 2.4), measured with integrated flux. We compare the star-forming main sequence with that from P. Popesso et al. (2023),<fnref n="28" /> who conducted a comprehensive literature search, and used a consistent method to derive the star-forming main sequence over $1 < z < 6$ for galaxies with $M_{*} > 10^{8.5} M_{\odot}$. As Figure 10 shows, our measurements are in excellent agreement with the distribution measured by P. Popesso et al. (2023). We also plot the star-forming main sequence from J. Leja et al. (2022), who also used PROSPECTOR with a very similar SED model to ours, but only measured

<footnote n="27">²⁷ We base this estimate on the stellar-to-halo mass ratio, and the halo mass growth history from P. S. Behroozi et al. (2013a, 2013b). For a $10^{9.7} M_{\odot}$ galaxy at $z = 3$, we find it is hosted by a $\sim 10^{12.7} M_{\odot}$ halo at $z \sim 0.1$, which is ≈0.7 dex more massive than the Milky Way’s dark matter halo (L. Posti & A. Helmi 2019).</footnote>

<footnote n="28">²⁸ As noted by P. Popesso et al. (2023), they made a –0.3 dex correction for the stellar-mass measures from the PROSPECTOR fitting with nonparametric SFHs. This systematic offset in stellar-mass measures between nonparametric and parametric SFHs has also been reported by others (J. Leja et al. 2019, 2022; Z. Ji & M. Giavalisco 2022, 2023). For a proper comparison, we therefore enlarge the stellar mass in the best-fit star-forming main sequence of P. Popesso et al. (2023; their Equation (10)) by 0.3 dex.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** The stellar-mass surface density at the location of C1. The symbol used for each galaxy is the same as in Figure 9. The solid line shows the Σ<sub>1 kpc</sub> evolution of galaxies at 0.5 < z < 3 measured by G. Barro et al. (2017). We extrapolate this evolution to z > 3 and plot it as the dashed line. The evolution of star-forming galaxies and quiescent galaxies is plotted in cyan and magenta, respectively. The shaded region marks the ±0.3 dex range. Although the six galaxies have not quenched yet, their stellar-mass densities within the region of C1 are already high, approaching the central stellar-mass surface density of high-redshift quiescent galaxies (see Section 4.2 for detailed discussions and caveats of this comparison).</figcaption></figure>

the relation for galaxies up to z = 2.7 with $M_* > 10^{10.2} M_\odot$ owing to the stellar-mass incompleteness. As Figure 10 shows, the six galaxies presented in this work are below the star-forming main sequence by 0.2–0.7 dex, suggesting a possible link between the presence of a massive stellar core (e.g., C1) and galaxy quenching.

In Figure 11, we continue the discussion by comparing the stellar-mass surface densities at the location of C1 ($\Sigma_{\text{C1}}$) with the prediction based on the evolution of $\Sigma_{\text{1kpc}}$, i.e., the stellar-mass surface density within the central radius of 1 kpc, of 0.5 < z < 3 galaxies measured by G. Barro et al. (2017) using HST data from CANDELS. Specifically, we compare with the $\Sigma_{\text{1 kpc}}$ evolution of star-forming (the cyan line in Figure 11) and quiescent (the magenta line in Figure 11) galaxies with $M_* = 10^{9.7} M_\odot$, i.e., the median stellar mass of the galaxies presented in this work.

Before proceeding to discuss the results, we clarify how we estimate $\Sigma_{\text{C1}}$ for a fair comparison with the existing $\Sigma_{\text{1 kpc}}$ measures, and caution about the caveats in this comparison. First, because C1 are close to the center of stellar-light distributions (NIRCam images at ~4µm in Figures 3 and 4), and their sizes are ≈0.5–1 kpc (effective radius, Figure 8), $\Sigma_{\text{C1}}$ is a good proxy for the central stellar-mass surface densities of the galaxies presented in this work. Second, we note that the existing $\Sigma_{\text{1 kpc}}$ measures simply assumed a single Sérsic profile for the entire galaxy, meaning that those $\Sigma_{\text{1 kpc}}$ measures included all stellar masses in the central regions, regardless of their origins (e.g., cores versus stars in disks). For a fair comparison, instead of using $M_{\text{C1}}^{\text{int}}$ (Equation (4)), we thus include all masses at the location of C1, i.e., $\Sigma_{\text{C1}} = \eta M_{\text{C1}}/A_{\text{C1}}$, where $A_{\text{C1}}$ is the aperture size used for C1. Third, the earlier studies were only able to study the evolution of $\Sigma_{\text{1 kpc}}$ up to z = 3. In this work, we therefore need to extrapolate the $\Sigma_{\text{1 kpc}}$ evolution from G. Barro et al. (2017) to a higher redshift range, i.e., 3 < z < 4.5. Finally, at z ~ 3, while the stellar-mass

completeness of CANDELS is still ~$10^{9.5} M_\odot$ for star-forming galaxies (e.g., Z. Ji et al. 2018; P. Santini et al. 2022), it rises to ~$10^{10.3} M_\odot$ for quiescent galaxies (G. Barro et al. 2017). We therefore also need to extrapolate the $\Sigma_{\text{1 kpc}}$ evolution for quiescent galaxies from G. Barro et al. (2017) to a lower stellar mass $10^{9.7} M_\odot$.

As Figure 11 shows, except JADES-JEMS-16296, all five other galaxies have large central stellar-mass surface densities, namely, that their $\Sigma_{\text{C1}}$ are already similar to the $\Sigma_{\text{1 kpc}}$ of quiescent galaxies of similar redshifts and masses. Interestingly, JADES-JEMS-16296 in fact is the galaxy that shows evidence of recent rejuvenation of central star formation, which will be further discussed in Section 4.3. Although the galaxies presented in this work are not quenched yet, they are below the star-forming main sequence. On the one hand, although we cannot be sure that these six galaxies will not rejuvenate their star formation activities in the future, the findings here suggest that galaxies have already formed significant amount of their mass in a spheroidal central component while they are still in the star-forming phase, which is consistent with what has been found in the IllustrisTNG simulations (A. Pillepich et al. 2019; S. Tacchella et al. 2019). On the other hand, it is also possible that the galaxies are actually heading toward quiescence. The presence of a massive stellar core with a surface density comparable to quenched galaxies then suggests that the suppression of star formation may occur along with the development of dense central regions. The implication is that there is either a causal link between quenching and central structural transformations, or an as-yet unidentified physical process that controls both events.

The conclusion above is in broad agreement with the studies at lower redshifts (e.g., J. J. Fang et al. 2013; G. Barro et al. 2017; K. E. Whitaker et al. 2017; S. Tacchella et al. 2018; K. A. Suess et al. 2021; Z. Ji & M. Giavalisco 2023). What is different here, though, is that we are now able to study the buildup of central regions in the act at 3 < z < 5 in the progenitor of lower-redshift massive quiescent galaxies, thanks to the deep NIRCam imaging obtained by JADES and JEMS. This, in our view, is an important step forward because those lower-redshift quiescent galaxies are already old—the typical stellar age of z ~ 1–2 massive quiescent galaxies is about 1 Gyr (e.g., A. C. Carnall et al. 2019b; S. Tacchella et al. 2022a; Z. Ji & M. Giavalisco 2022). This fact complicates the interpretation for the correlation between the star formation properties and central stellar-mass surface density that is well established at z = 1–2 (e.g., S. Belli et al. 2015; C. C. Williams et al. 2017; V. Estrada-Carpenter et al. 2020), because post-quenching processes can alter galaxies’ structures. Thus, it is difficult to assess if the correlation is really due to a causal link between the development of dense cores and quenching, or is simply a byproduct of another yet-to-be identified physical process. Pushing observations closer to higher redshifts z > 3, i.e., closer to the time of quenching, enables a more direct view of the relative timing between structural transformations and quenching, hence, guarantees a more accurate constraint on the relationship between the two phenomena.

### 4.3. Consistency between Our Observations and the Gas-rich Compaction Scenario

Finally, we discuss our constraints on the physical mechanisms responsible for the correlation between the star formation properties and central stellar-mass surface density observed in
galaxies at cosmic noon (e.g., M. Franx et al. 2008; E. Cheung et al. 2012; G. Barro et al. 2017; K. E. Whitaker et al. 2017). Our observations are fully in line with the so-called “gas-rich compaction” scenario, which, as we already mentioned in Section 1, is a class of processes whose main feature is highly dissipative gas accretion toward the center of galaxies. This process promotes enhanced central star formation activity, followed by the central gas depletion and finally the quenching of central star formation (A. Dekel et al. 2009; D. Ceverino et al. 2010; A. Dekel & A. Burkert 2014; A. Zolotov et al. 2015; S. Wellons et al. 2015; S. Tacchella et al. 2016).

One important prediction of gas-rich compaction is that, right after a compaction event, we should observe a younger stellar population at the center as a result of recent central starburst triggered by the highly dissipative gas accretion. This is observed in the six galaxies where C1 is comparable to/younger than the smooth component (Section 3.1). Additionally, we find that C1 features a common SFH whose SFR declines over the past 0.1–0.3 Gyr, following a major star-forming episode ≈0.5–1 Gyr prior, and peaked ≈0.2–0.5 Gyr ago. This SFH is consistent with the six galaxies experiencing a gas compaction event, followed by a suppression of star formation likely as the result of central gas depletion. If there is future gas replenishment, this chain of events, i.e., gas compaction followed by gas depletion, will repeat several times before the final quiescence of the entire galaxy at $z < 3$ (S. Tacchella et al. 2016). Indeed, we see evidence of star formation rejuvenation in the C1 of JADES-JEMS-16296 by the time of observation. This rejuvenation event happens at ≈0.3 Gyr after the peak of the galaxy’s most recent major star formation episode. Interestingly, this timescale is very similar to the model predictions for the oscillation timescale (0.2–0.5 $t_H$, corresponding to 0.3–0.8 Gyr at $z=3.7$) of star formation triggered by gas replenishment (S. Tacchella et al. 2016). However, we caution that this result is sensitive to the assumed prior of nonparametric SFHs (see Appendix C).

## 5. Summary

This is the first paper of the series that uses deep 14-filter JWST/NIRCam imaging data from the JEMS and JADES surveys to study in detail the spatially resolved stellar populations of galaxies at redshifts $3 < z < 4.5$, a key epoch of the Universe that can provide unique insight into the role that substructures among massive galaxies may play in the transition from star-forming disks to bulge-dominated, quiescent galaxies.

We visually inspected color maps at rest-frame 3600–4100 Å produced by JWST/NIRCam imaging for all $3 < z < 4.5$ galaxies in JEMS, from which we identified a sample of 37 galaxies having substructures with distinct rest-frame colors. These 37 galaxies show diverse spatially resolved stellar populations. In this first paper, we present the spatially resolved stellar populations in six galaxies with massive stellar cores, representing a specific stage of galaxy structural transformation.

We divided each one of the six galaxies into different structural components, including (1) a bright central stellar core C1, the focus of this study, (2) off-center, minor clumps C2 and C3 (if present), and (3) the smooth component representing the remaining, extended light distribution of the galaxy. Together with nine-filter ancillary imaging from HST, we fit a total of

23-filter photometry covering truly panchromatic swathes of observed wavelength range $\lambda = 0.4–5\mu$m, corresponding to the rest-frame 0.1–1$\mu$m, using the SED fitting code PROSPECTOR with nonparametric SFHs. We then studied the stellar population properties of those different structural components.

In each one of the galaxies, C1 has a redder color at rest-frame 4000 Å and a significantly lower sSFR than off-center clumps C2 and C3. We showed that C1 is $\gtrsim 2$ times more massive than $M_{\text{Toomre}}$, i.e., the maximum stellar mass that a clump can form via Toomre instability, suggesting that it may not have formed via single in situ fragmentation in an unstable, gaseous disk. Instead, because C1 has a typical stellar age of ~0.5 Gyr, which is similar to the timescale of clump inward migration as the result of dynamical friction, we argued that C1 likely formed through the coalescence of giant stellar clumps, a model prediction from the VDI, which is one of the processes associated with the gas-rich compaction scenario at high redshift.

We found that all six galaxies are below the star-forming main sequence by 0.2–0.7 dex. We showed that the stellar-mass surface densities of C1 are already comparable to the central stellar-mass surface density of quenched galaxies of similar masses and redshifts. In addition, we found that the stellar populations of C1 are either comparable to or younger than the smooth component. Putting these together, we thus concluded that we are likely witnessing the coeval development of dense central cores, along with the beginning of galaxy-wide quenching at $z > 3$, likely triggered by gas-rich compaction.

As a closing remark, we stress that the sample presented here only contains six galaxies and hence has no statistical power. In addition, the visual selection means that this study is biased toward bright substructures in relatively bright galaxies. To obtain a robust statistical characterization of massive stellar cores in high-redshift galaxies, and a comprehensive picture of its relationship with quenching, a much larger sample with a more uniform selection method is required. We, however, view this work as a pilot study to demonstrate the power of JWST, especially its imaging capability with NIRCam, in constraining the spatially resolved mass assembly history of galaxies at $z > 3$ in unprecedented detail. We reiterate that the galaxies here were selected from JEMS, whose sky coverage is relatively small (~10 arcmin²). The ongoing JADES survey will eventually produce NIRCam imaging and NIRSpec spectroscopy of similar depth over a sky area of $\gtrsim 200$ arcmin² in the GOODS fields, promising a statistically significant galaxy sample at $z > 3$ that will dramatically advance our understanding of the physics of galaxy structural transformation and quenching in the early Universe.

## Acknowledgments

Z.J., B.D.J., B.R., F.S., D.J.E., M.R., G.R., K.H., C.N.A.W., E.E., Z.C., J.M.H., and J.L. acknowledge funding from JWST/NIRCam contract to the University of Arizona NAS5-02015. The research of C.C.W. is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. W.M.B., T.J.L., R.M., and L.S. acknowledge support by the Science and Technology Facilities Council (STFC), ERC Advanced grant 695671 “QUENCH.” A.J.B. and J.C. acknowledge funding from the “FirstGalaxies” Advanced Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (grant agreement No. 789056). D.J.E. is supported as a Simons Investigator. R.H. acknowledges funding for this research was provided by the Johns Hopkins University, Institute for Data Intensive Engineering and Science (IDIES). S.C. acknowledges support by European Union's HE ERC Starting grant No. 101040227—WINGS. E. C.L. acknowledges support of an STFC Webb Fellowship (ST/ W001438/1). R.M. acknowledges funding from a research professorship from the Royal Society. This research of K.B. is supported in part by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project No. CE170100013.

This work is based on observations, JADES (M. Rieke et al. 2023a) and JEMS (C. Williams et al. 2023a), made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. This research is based (in part) on observations made with the NASA/ESA HST obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

This material is based upon High Performance Computing (HPC) resources supported by the University of Arizona TRIF, UITS, and Research, Innovation, and Impact (RII) and maintained by the University of Arizona Research Technologies department. This work made use of the lux supercomputer at UC Santa Cruz, which is funded by NSF MRI grant AST 1828315, as well as the High Performance Computing (HPC) resources at the University of Arizona, which is funded by the Office of Research Discovery and Innovation (ORDI), Chief Information Officer (CIO), and University Information Technology Services (UITS).

Facilities: HST, JWST.

*Software:* Prospector (B. D. Johnson et al. 2021), FSPS (C. Conroy et al. 2009; C. Conroy & J. E. Gunn 2010), MIST (A. Dotter 2016; J. Choi et al. 2016), MILES (J. Falcón-Barroso et al. 2011), Webbpsf (M. D. Perrin et al. 2012, 2014), Photutils (L. Bradley et al. 2022).

# Appendix A

## Comparing Effective PSF with Webbpsf-Predicted Model PSF

Knowledge of the point-spread function (PSF) plays an essential role in getting accurate PSF-matched aperture photometry and, hence, the inferences of physical properties of stellar populations. Per JWST User Documentation,<fnref n="29" /> NIRCam’s PSFs are under-sampled in both <2µm SW filters and <4µm LW filters. To overcome this, we have constructed empirical and theoretical PSF models for JADES and JEMS with the following three different approaches:

1. Effective PSF (ePSF)—We construct ePSFs using the EPSFBuilder in PHOTUTILS, which adopts the empirical method from J. Anderson & I. R. King (2000) to iteratively solve the centroids and fluxes of a list of input point sources, and then stack them together. For the list of point sources, we visually selected 13 isolated point sources from JADES and JEMS. These point sources are bright, but unsaturated across all NIRCam filters used in this study.

2. Model PSF (mPSF)—Two effects are taken into account during the construction of mPSFs. The first effect comes from instrumental responses, which are taken care by the software WEBBPSF. The other effect comes from data reduction, especially from the mosaicking process. To include both effects to the mPSF construction, we assume a list of synthetic point sources with known WCS coordinates. For each one of the point sources, we first project its WCS coordinate to individual exposures. In each exposure, we then convert the WCS coordinate to the corresponding detector coordinate, and use the WEBBPSF to predict the PSF at that location, during which we adopt the JWST on-orbit measured optical path difference map closest to the date of that exposure. We repeat these for all synthetic point sources. Finally, we use the mosaicking pipeline of JADES to combine individual exposures together and get the mosaicked PSF images. We have made the PSF mosaics for all NIRCam filters used in this work (three examples are shown in Figure 12). These mosaics allow us to check astrometric accuracy, image quality, and to study the spatial variation of PSFs, which will be presented in detail in Z. Ji et al. (2023, in preparation).

3. Simple WEBBPSF—Unlike mPSFs, simple WEBBPSF models do not include any observational and data reduction effects.

In Figures 13 and 14, we compare the above three PSF models. Excellent agreement is seen between ePSFs and mPSFs, with a typical difference of ≲1% in radial profiles of enclosed energy. We have tested our PSF-matched photometry using the mPSFs, and found no substantial changes at all for the results of this paper.

We note, however, that the analysis here demonstrates that simple WEBBPSF models can significantly deviate from ePSFs, especially at small angular scales of <0.′′1. Therefore, using simple WEBBPSF models can lead to significant systematics in the PSF-related analysis, if observational and data reduction effects are not taken into account. Once these effects are properly modeled, the prediction of WEBBPSF is very accurate. This last point is very important for extragalactic deep surveys like JADES, whose observational configurations (e.g., dithering patterns) can be complex and where not enough high-S/N, unsaturated stars are in the field for constructing high-quality ePSF models to large angular scales.

<footnote n="29">²⁹ https://jwst-docs.stsci.edu</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** JADES PSF mosaics of NIRCam/F090W (left), F200W (middle), and F444W (right) filters produced following the procedure described in Appendix A.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Comparisons among three different PSF models (Appendix A) for the NIRCam filters used by JADES. Plotted in the bottom of each panel are the differences between mPSFs and ePSFs (red dashed line), and between simple WEBBPSF models and ePSFs (blue dashed line), where the horizontal, gray dashed lines mark ±1% range. Excellent agreement is observed between ePSFs and mPSFs, while simple WEBBPSF models significantly deviate from ePSFs at small angular scales.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14. Similar to Figure 13, but for the five medium-band NIRCam filters used by JEMS.</figcaption></figure>

## Appendix B

### Comparing Coadded SFHs with SFHs from SED Fitting of Integrated Flux

Apart from comparing the coadded stellar mass and SFR with those measured from SED fitting of integrated flux (Figure 5), here we compare the coadded SFH of individual

structural components, with the SFH from SED fitting of integrated flux. The results are shown in Figure 15. The two SFHs are consistent with each other within uncertainty, further demonstrating the consistency and robustness of our SED fitting results.

<figure><img src="image_16.png" /><figcaption>Figure 15. Comparison of the coadded SFH of individual structural components with the SFH derived from SED fitting of integrated flux. The y-axis shows the relative difference, i.e., (coadded—integrated) / integrated, of the two SFHs. The gray shaded region shows $1\sigma$ uncertainty.</figcaption></figure>
## Appendix C
### SED Fitting with Different Priors of Star Formation History

Apart from the fiducial SED fitting with nonparametric SFHs assuming the continuity prior (Section 2.4), we have also tested our PROSPECTOR fitting results using the other two SFHs, namely, nonparametric SFHs with the Dirichlet prior (J. Leja et al. 2017) and parametric SFHs with a delayed-tau functional form. For the difference between the continuity and Dirichlet priors, we refer readers to J. Leja et al. (2017) for details. In short, these two priors for nonparametric SFHs are very similar except that the continuity prior is strongly against sudden changes in SFR in adjacent lookback time bins. For the delayed-tau model, i.e., $SFR(t) \propto (t - t_{age})e^{-(t-t_{age})/\tau}$, we fit $\tau$ with a logarithmically

flat prior of $log(\tau/Gyr) \in (-2, 1)$, and $t_{age}$ with a flat prior of $t_{age}/Gyr \in (10^{-3}, t_{H})$.

In Figure 16, we compare the measurements of stellar mass and SFR from different SFH assumptions. Very good agreement is seen across the three SFH models. In Figures 17 and 18, we compare the fiducial reconstructed SFHs with the delayed-tau and Dirichlet nonparametric SFHs, respectively. As expected, the nonparametric SFHs are able to capture more complex shapes than the delayed-tau ones, but the overall shapes of these two SFHs are very similar. Meanwhile, relative to the fiducial SFHs from the continuity prior, the Dirichlet nonparametric SFHs show larger fluctuations in SFR in adjacent lookback time bins, but these two nonparametric SFHs, in general, mirror each other in terms of their overall shapes. The tests here show that the results of this work are not sensitive to the assumptions of SED models.

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Comparisons of stellar-mass and SFR measures from different SFHs for the individual structural components of the six galaxies presented in this work. Color-coding is the same as Figure 6 in the main text (C1 are in red). Symbols used for individual galaxies are the same as Figure 9 in the main text. The black solid line in each panel marks the one-to-one relationship.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 17.** The comparison of the nonparametric SFH reconstructed by our fiducial model (Section 2.4), plotted as dashed lines, with the SFH derived from a parametric delayed-tau model, plotted as solid lines (the best-fit) and shaded regions ($1\sigma$ uncertainty).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Similar to Figure 17, but now comparing with the nonparametric SFH reconstructed using the Dirichlet prior.</figcaption></figure>

## Appendix D
Detailed Descriptions of the Spatially Resolved Stellar Populations of the Six Galaxies

JADES-JEMS-13396. The top row of Figure 6 shows the SED fitting result of this galaxy. Despite that the mass-weighted stellar age of C1 is consistent with that of the smooth

component, the shapes of their SFHs are different. For the smooth component, its mass assembly was at a relatively low, but nonnegligible, rate in early times of $t_{\text{lookback}} > 0.4$ Gyr, followed by a rapid increase in SFR. For C1, however, most of its mass formed through a recent star formation episode during $t_{\text{lookback}} = 0.2 \sim 0.8$ Gyr. By the time of observation, the smooth component is forming stars at a higher rate than C1,
with the sSFR of the former being ≈1.6 times larger (Table 2). These differences in SFH shapes qualitatively remain if we instead use a parametric delayed-tau or nonparametric Dirichlet SFH priors (Figures 17 and 18). We also notice that a color gradient is clearly present in the smooth component, with the upper-left side (likely contains several faint clumps in it) being bluer than its rest parts.

JADES-JEMS-15157. The middle row of Figure 6 shows the SED fitting result of this galaxy. Unlike C2 having a monotonically rising SFH, the SFH of C1 gradually rises toward a relatively recent peak at $t_{\text{lookback}} \approx 0.3$ Gyr, followed by a rapid decline over the past 0.2 Gyr. By the time of observation, the sSFR of C1 is only $\approx 0.4 \text{ Gyr}^{-1}$, $\approx 7.5$ times smaller than that of C2 (Table 2). We note that C3 has a similar stellar population to that of C1, but this is because the spatial separation between C1 and C3 is so small—about 0.′′1, which is smaller than the angular resolution of NIRCam/F444W imaging—that the aperture photometry (Section 2.3) for one substructure is significantly affected by the other. Finally, while C1 and the smooth component share a similar SFH in early times, they differentiate from each other at late times: unlike C1, whose SFR quickly declined after the $t_{\text{lookback}} \approx 0.3$ Gyr peak, the smooth component retained that peak SFR toward a later time of $t_{\text{lookback}} < 0.1$ Gyr. All of these remain qualitatively unchanged if we instead use a parametric delayed-tau or nonparametric Dirichlet SFH priors (Figures 17 and 18).

JADES-JEMS-6885. The bottom row of Figure 6 shows the SED fitting result of this galaxy. There is no significant mass assembly before $t_{\text{lookback}} > 1$ Gyr for all structural components, but their SFHs are different over the recent 1 Gyr. Compared to C2 and C3, C1 is older by 0.3–0.5 Gyr (Table 2). Moreover, C1 has a post-starburst-like SFH; namely, it experienced a burst of star formation during $t_{\text{lookback}} = 0.2–0.6$ Gyr, followed by a rapid decline of SFR. Differing from C1, both C2 and C3 are experiencing ongoing bursts of star formation that started $\approx 0.2$ Gyr ago. By the time of observation, the sSFR of C1 is 6.3 and 9.3 times smaller than those of C2 and C3, respectively. Finally, while the smooth component generally mirrored the shape of C1’s SFH before the peak ($t_{\text{lookback}} \approx 0.2$ Gyr), they differed from each other afterward: the SFR of C1 quickly declined to a very low level, while the smooth component retained the peak SFR to a later time of $t_{\text{lookback}} < 0.1$ Gyr. All of these differences in the SFHs of different structural components remain qualitatively unchanged if we instead use the other two SFH priors (Figures 17 and 18).

JADES-JEMS-16296. The top row of Figure 7 shows the SED fitting result of this galaxy. Both C1 and C2 are experiencing bursts of star formation started $\lesssim 0.1$ Gyr ago, and their instantaneous sSFRs are similar, $\approx 5 \text{ Gyr}^{-1}$. While it is not seen from the delayed-tau SFH (Figure 17), the reconstructed nonparametric SFH (Figures 7 and 18) suggests that older ($t_{\text{lookback}} \gtrsim 0.5$ Gyr) stellar populations are also present in C1, leading to its larger mass-weighted stellar age than C2 (Table 2).

JADES-JEMS-14436. The middle row of Figure 7 shows the SED fitting result of this galaxy. The major mass assembly event of C1 started at $t_{\text{lookback}} \approx 0.6$ Gyr, peaked at $t_{\text{lookback}} \approx 0.2$ Gyr, and gradually declined afterward. For C2, its SFR monotonically increased since $t_{\text{lookback}} \approx 0.8$ Gyr, before which its SFR was consistent with zero. There also seems to be a tidal tail associated with C2 (see, e.g., NIRCam F150W cutouts in Figure 4), so this is a possible merging system. Also noticed for

this galaxy is that the smooth component has a very extended SFH, namely, its mass assembly began at $t_{\text{lookback}} > 1$ Gyr, corresponding to redshift $z > 7$, and the SFR started to decline since $t_{\text{lookback}} \approx 0.6$ Gyr. All of these features in SFHs remain qualitatively unchanged if we instead use the other two SFH priors (Figures 17 and 18).

JADES-JEMS-11059. The bottom row of Figure 7 shows the SED fitting result of this galaxy. C1 contains old stellar populations with age $> 1$ Gyr, and it also has a recent burst of star formation that started $\approx 0.2$ Gyr ago, peaked at $t_{\text{lookback}} \approx 0.1$ Gyr, and declined afterward. Meanwhile, C3 has an overall rising SFH. For C2, although it is not in the segmentation map of the galaxy (Figure 4), it is still included to the analysis because its photometric redshift solution is consistent with that of the galaxy JADES-JEMS-11059. C2 is experiencing a starburst started $\lesssim 0.3$ Gyr ago.

## ORCID iDs

Zhiyuan Ji https://orcid.org/0000-0001-7673-2257
Christina C. Williams https://orcid.org/0000-0003-2919-7495

Sandro Tacchella https://orcid.org/0000-0002-8224-4505

Katherine A. Suess https://orcid.org/0000-0002-1714-1905

William M. Baker https://orcid.org/0000-0003-0215-1104

Stacey Alberts https://orcid.org/0000-0002-8909-8782

Andrew J. Bunker https://orcid.org/0000-0002-8651-9879

Benjamin D. Johnson https://orcid.org/0000-0002-9280-7594

Brant Robertson https://orcid.org/0000-0002-4271-0364

Fengwu Sun https://orcid.org/0000-0002-4622-6617

Daniel J. Eisenstein https://orcid.org/0000-0002-2929-3121

Marcia Rieke https://orcid.org/0000-0002-7893-6170

Michael V. Maseda https://orcid.org/0000-0003-0695-4414

Kevin Hainline https://orcid.org/0000-0003-4565-8239

Ryan Hausen https://orcid.org/0000-0002-8543-761X

George Rieke https://orcid.org/0000-0003-2303-6519

Christopher N. A. Willmer https://orcid.org/0000-0001-9262-9997

Eiichi Egami https://orcid.org/0000-0003-1344-9475

Irene Shivaei https://orcid.org/0000-0003-4702-7561

Stefano Carniani https://orcid.org/0000-0002-6719-380X

Stephane Charlot https://orcid.org/0000-0003-3458-2275

Jacopo Chevallard https://orcid.org/0000-0002-7636-0534

Emma Curtis-Lake https://orcid.org/0000-0002-9551-0534

Tobias J. Looser https://orcid.org/0000-0002-3642-2446

Roberto Maiolino https://orcid.org/0000-0002-4985-3819

Chris Willott https://orcid.org/0000-0002-4201-7367

Zuyi Chen https://orcid.org/0000-0002-2178-5471

Jakob M. Helton https://orcid.org/0000-0003-4337-6211

Jianwei Lyu https://orcid.org/0000-0002-6221-1829

Erica Nelson https://orcid.org/0000-0002-7524-374X

Rachana Bhatawdekar https://orcid.org/0000-0003-0883-2226

Kristan Boyett https://orcid.org/0000-0003-4109-304X

Lester Sandles https://orcid.org/0000-0001-9276-7062

## References

* Anderson, J., & King, I. R. 2000, PASP, 112, 1360

* Athanassoula, E., Laurikainen, E., Salo, H., & Bosma, A. 2015, MNRAS, 454, 3843

* Bacon, R., Brinchmann, J., Conseil, S., et al. 2022, A&A, 670, 42
Balogh, M. L., Morris, S. L., Yee, H. K. C., Carlberg, R. G., & Ellingson, E. 1999, ApJ, 527, 54

Barbuy, B., Chiappini, C., & Gerhard, O. 2018, ARA&A, 56, 223

Barnes, J. E., & Hernquist, L. 1996, ApJ, 471, 115

Barro, G., Faber, S. M., Koo, D. C., et al. 2017, ApJ, 840, 47

Behrendt, M., Burkert, A., & Schartmann, M. 2015, MNRAS, 448, 1007

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013a, ApJ, 770, 57

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013b, ApJL, 762, L31

Belli, S., Newman, A. B., & Ellis, R. S. 2015, ApJ, 799, 206

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton, NJ: Princeton Univ. Press)

Boucaud, A., Bocchio, M., Abergel, A., et al. 2016, A&A, 596, A63

Bournaud, F. 2016, in Galactic Bulges, ed. E. Laurikainen, R. Peletier, & D. Gadotti (Berlin: Springer), 355

Bournaud, F., Elmegreen, B. G., & Elmegreen, D. M. 2007, ApJ, 670, 237

Bournaud, F., Jog, C. J., & Combes, F. 2005, A&A, 437, 69

Bournaud, F., Perret, V., Renaud, F., et al. 2014, ApJ, 780, 57

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

Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151

Bruce, V. A., Dunlop, J. S., Cirasuolo, M., et al. 2012, MNRAS, 427, 1666

Bruce, V. A., Dunlop, J. S., McLure, R. J., et al. 2014, MNRAS, 444, 1001

Bruzual, A. G. 1983, ApJ, 273, 105

Buck, T., Macciò, A. V., Obreja, A., et al. 2017, MNRAS, 468, 3628

Bundy, K., Scarlata, C., Carollo, C. M., et al. 2010, ApJ, 719, 1969

Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2022, JWST Calibration Pipeline, v1.8.1, Zenodo, doi:10.5281/ZENODO.7215690

Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019a, ApJ, 873, 44

Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2019b, MNRAS, 490, 417

Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2023a, Natur, 7971, 716

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

Catinella, B., Schiminovich, D., Cortese, L., et al. 2013, MNRAS, 436, 34

Ceverino, D., Dekel, A., & Bournaud, F. 2010, MNRAS, 404, 2151

Ceverino, D., Dekel, A., Mandelker, N., et al. 2012, MNRAS, 420, 3490

Ceverino, D., Klypin, A., Klimek, E. S., et al. 2014, MNRAS, 442, 1545

Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718

Cheung, E., Faber, S. M., Koo, D. C., et al. 2012, ApJ, 760, 131

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833

Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486

Conselice, C. J. 2014, ARA&A, 52, 291

Conselice, C. J., Bershady, M. A., & Jangren, A. 2000, ApJ, 529, 886

Cowie, L. L., Hu, E. M., & Songaila, A. 1995, AJ, 110, 1576

Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686

Danovich, M., Dekel, A., Hahn, O., Ceverino, D., & Primack, J. 2015, MNRAS, 449, 2087

Debattista, V. P., Carollo, C. M., Mayer, L., & Moore, B. 2004, ApJL, 604, L92

Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870

Dekel, A., & Krumholz, M. R. 2013, MNRAS, 432, 455

Dekel, A., Sari, R., & Ceverino, D. 2009, ApJ, 703, 785

Dessauges-Zavadsky, M., & Adamo, A. 2018, MNRAS, 479, L118

Dimauro, P., Daddi, E., Shankar, F., et al. 2022, MNRAS, 513, 256

Dotter, A. 2016, ApJS, 222, 8

Eisenstein, D. J., Ferruit, P., Rieke, M. J., Willmer, C. N. A., & Willott, C. J. 2017, arXiv:2306.02465

Eliche-Moral, M. C., González-García, A. C., Balcells, M., et al. 2011, A&A, 533, A104

Elmegreen, B. G., Bournaud, F., & Elmegreen, D. M. 2008, ApJ, 688, 67

Elmegreen, B. G., & Elmegreen, D. M. 2005, ApJ, 627, 632

Elmegreen, D. M., Elmegreen, B. G., Marcus, M. T., et al. 2009, ApJ, 701, 306

Elmegreen, D. M., Elmegreen, B. G., Ravindranath, S., & Coe, D. A. 2007, ApJ, 658, 763

Escala, A., & Larson, R. B. 2008, ApJL, 685, L31

Estrada-Carpenter, V., Papovich, C., Momcheva, I., et al. 2020, ApJ, 898, 171

Falcón-Barroso, J., Sánchez-Blázquez, P., Vazdekis, A., et al. 2011, A&A, 532, A95

Fang, J. J., Faber, S. M., Koo, D. C., & Dekel, A. 2013, ApJ, 776, 63

Ferguson, H. C., Dickinson, M., Giavalisco, M., et al. 2004, ApJL, 600, L107

Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, ApJ, 706, 1364

Förster Schreiber, N. M., Genzel, R., Lehnert, M. D., et al. 2006, ApJ, 645, 1062

Förster Schreiber, N. M., Renzini, A., Mancini, C., et al. 2018, ApJS, 238, 21

Förster Schreiber, N. M., Shapley, A. E., Genzel, R., et al. 2011, ApJ, 739, 45

Förster Schreiber, N. M., & Wuyts, S. 2020, ARA&A, 58, 661

Franx, M., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2008, ApJ, 688, 770

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1

Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 20

Genel, S., Naab, T., Genzel, R., et al. 2012, ApJ, 745, 11

Genzel, R., Burkert, A., Bouché, N., et al. 2008, ApJ, 687, 59

Genzel, R., Förster Schreiber, N. M., Lang, P., et al. 2014, ApJ, 785, 75

Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101

Genzel, R., Price, S. H., Übler, H., et al. 2020, ApJ, 902, 98

Genzel, R., Tacconi, L. J., Eisenhauer, F., et al. 2006, Natur, 442, 786

Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20

Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJL, 600, L93

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Guedes, J., Mayer, L., Carollo, M., & Madau, P. 2013, ApJ, 772, 36

Guo, Y., Ferguson, H. C., Bell, E. F., et al. 2015, ApJ, 800, 39

Guo, Y., Giavalisco, M., Ferguson, H. C., Cassata, P., & Koekemoer, A. M. 2012, ApJ, 757, 120

Guo, Y., Rafelski, M., Bell, E. F., et al. 2018, ApJ, 853, 108

Hasselquist, S., Zasowski, G., Feuillet, D. K., et al. 2020, ApJ, 901, 109

Hernquist, L. 1989, Natur, 340, 687

Hopkins, P. F., Cox, T. J., Younger, J. D., & Hernquist, L. 2009, ApJ, 691, 1168

Hopkins, P. F., Kereš, D., Murray, N., Quataert, E., & Hernquist, L. 2012, MNRAS, 427, 968

Huertas-Company, M., Gravet, R., Cabrera-Vives, G., et al. 2015, ApJS, 221, 8

Huertas-Company, M., Guo, Y., Ginzburg, O., et al. 2020, MNRAS, 499, 814

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Ilbert, O., Salvato, M., Le Floc'h, E., et al. 2010, ApJ, 709, 644

Illingworth, G., Magee, D., Bouwens, R., et al. 2016, arXiv:1606.00841

Immeli, A., Samland, M., Westera, P., & Gerhard, O. 2004, ApJ, 611, 20

Inoue, S., Dekel, A., Mandelker, N., et al. 2016, MNRAS, 456, 2052

Inoue, S., & Saitoh, T. R. 2012, MNRAS, 422, 1902

Ji, Z., & Giavalisco, M. 2022, ApJ, 935, 120

Ji, Z., & Giavalisco, M. 2023, ApJ, 943, 54

Ji, Z., Giavalisco, M., Kirkpatrick, A., et al. 2022, ApJ, 925, 74

Ji, Z., Giavalisco, M., Williams, C. C., et al. 2018, ApJ, 862, 135

Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22

Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2012, ApJ, 758, 106

Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 33

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Kormendy, J., & Kennicutt, R. C. J. 2004, ARA&A, 42, 603

Kriek, M., & Conroy, C. 2013, ApJL, 775, L16

Kroupa, P. 2001, MNRAS, 322, 231

Labbé, I., Huang, J., Franx, M., et al. 2005, ApJL, 624, L81

Lang, P., Wuyts, S., Somerville, R. S., et al. 2014, AJ, 788, 11

Law, D. R., Steidel, C. C., Erb, D. K., et al. 2009, ApJ, 697, 2057

Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3

Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170

Leja, J., Speagle, J. S., Ting, Y.-S., et al. 2022, ApJ, 936, 165

Liu, D., Schinnerer, E., Groves, B., et al. 2019, ApJ, 887, 235

Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97

Lotz, J. M., Primack, J., & Madau, P. 2004, AJ, 128, 163

Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2

Lyu, J., Alberts, S., Rieke, G. H., & Rujopakarn, W. 2022, ApJ, 941, 29

Madau, P. 1995, ApJ, 441, 18

Mandelker, N., Dekel, A., Ceverino, D., et al. 2014, MNRAS, 443, 3675

Mandelker, N., Dekel, A., Ceverino, D., et al. 2017, MNRAS, 464, 635

McGaugh, S. S., & Schombert, J. M. 2014, AJ, 148, 77

Merlin, E., Castellano, M., Santini, P., et al. 2021, A&A, 649, A22

Merlin, E., Fortuni, F., Torelli, M., et al. 2019, MNRAS, 490, 3309

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

Murata, K. L., Kajisawa, M., Taniguchi, Y., et al. 2014, ApJ, 786, 15

Murray, N., Quataert, E., & Thompson, T. A. 2010, ApJ, 709, 191

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18

Nelson, E., van Dokkum, P., Franx, M., et al. 2014, Natur, 513, 394

Nenkova, M., Sirocky, M. M., Ivezić, Ž., & Elitzur, M. 2008a, ApJ, 685, 147

Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezić, Ž., & Elitzur, M. 2008b, ApJ, 685, 160

Newman, A. B., Belli, S., Ellis, R. S., & Patel, S. G. 2018, ApJ, 862, 126

Newman, A. B., Belli, S., & Ellis, R. S. 2015, ApJL, 813, L7

Noguchi, M. 1999, ApJ, 514, 77
Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793

Oesch, P. A., Brammer, G., Naidu, R. P., et al. 2023, MNRAS, 525, 2864

Okamoto, T. 2013, MNRAS, 428, 718

Perez, J., Valenzuela, O., Tissera, P. B., & Michel-Dansac, L. 2013, MNRAS, 436, 259

Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363

Perret, V., Renaud, F., Epinat, B., et al. 2014, A&A, 562, A1

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

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

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6

Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526

Posti, L., & Helmi, A. 2019, A&A, 621, A56

Prochaska, J. X., & Wolfe, A. M. 2009, ApJ, 696, 1543

Quadri, R., Marchesini, D., van Dokkum, P., et al. 2007, AJ, 134, 1103

Rao, S. M., Turnshek, D. A., & Nestor, D. B. 2006, ApJ, 636, 610

Rieke, M., Robertson, B., Tacchella, S., et al. 2023a, ApJS, 269, 16

Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023b, PASP, 135, 028001

Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2022, NatAs, 7, 611

Saintonge, A., Kauffmann, G., Kramer, C., et al. 2011, MNRAS, 415, 32

Santini, P., Castellano, M., Fontana, A., et al. 2022, ApJ, 940, 135

Sattari, Z., Mobasher, B., Chartab, N., et al. 2023, ApJ, 951, 13

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

Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150

Shahidi, A., Mobasher, B., Nayyeri, H., et al. 2020, ApJ, 897, 44

Shapiro, K. L., Genzel, R., Förster Schreiber, N. M., et al. 2008, ApJ, 682, 231

Shibuya, T., Ouchi, M., & Harikane, Y. 2015, ApJS, 219, 15

Shibuya, T., Ouchi, M., Kubo, M., & Harikane, Y. 2016, ApJ, 821, 72

Simons, R. C., Kassin, S. A., Weiner, B. J., et al. 2017, ApJ, 843, 46

Skilling, J. 2004, in AIP Conf. Ser. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint (Melville, NY: AIP), 395

Soto, E., de Mello, D. F., Rafelski, M., et al. 2017, ApJ, 837, 6

Speagle, J. S. 2020, MNRAS, 493, 3132

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

Suess, K. A., Leja, J., Johnson, B. D., et al. 2022, ApJ, 935, 146

Tacchella, S., Carollo, C. M., Förster Schreiber, N. M., et al. 2018, ApJ, 859, 56

Tacchella, S., Carollo, C. M., Renzini, A., et al. 2015, Sci, 348, 314

Tacchella, S., Conroy, C., Faber, S. M., et al. 2022a, ApJ, 926, 134

Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 458, 242

Tacchella, S., Diemer, B., Hernquist, L., et al. 2019, MNRAS, 487, 5416

Tacchella, S., Eisenstein, D. J., Hainline, K., et al. 2023, ApJ, 952, 15

Tacchella, S., Finkelstein, S. L., Bagley, M., et al. 2022b, ApJ, 927, 170

Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Natur, 463, 781

Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179

Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246

Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157

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

Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2014, ApJ, 783, 85

Toomre, A. 1964, ApJ, 139, 1217

van den Bergh, S., Abraham, R. G., Ellis, R. S., et al. 1996, AJ, 112, 359

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

van der Wel, A., Rix, H.-W., Wuyts, S., et al. 2011, ApJ, 730, 38

Walter, F., Carilli, C., Neeleman, M., et al. 2020, ApJ, 902, 111

Weiner, B. J., Willmer, C. N. A., Faber, S. M., et al. 2006, ApJ, 653, 1027

Wellons, S., Torrey, P., Ma, C.-P., et al. 2015, MNRAS, 449, 361

Whitaker, K. E., Ashas, M., Illingworth, G., et al. 2019, ApJS, 244, 16

Whitaker, K. E., Bezanson, R., van Dokkum, P. G., et al. 2017, ApJ, 838, 19

Whitaker, K. E., Labbé, I., van Dokkum, P. G., et al. 2011, ApJ, 735, 86

Williams, C., Tacchella, S., & Maseda, M. 2023a, Data from the JWST Extragalactic Medium-band Survey (JEMS), STScI/MAST

Williams, C. C., Giavalisco, M., Bezanson, R., et al. 2017, ApJ, 838, 94

Williams, C. C., Giavalisco, M., Cassata, P., et al. 2014, ApJ, 780, 1

Williams, C. C., Tacchella, S., Maseda, M. V., et al. 2023b, ApJS, 268, 15

Wisnioski, E., Förster Schreiber, N. M., Fossati, M., et al. 2019, ApJ, 886, 124

Wuyts, S., Förster Schreiber, N. M., Genzel, R., et al. 2012, ApJ, 753, 114

Xiang, M., & Rix, H.-W. 2022, Natur, 603, 599

Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327