The MAGE Project: Size Evolution and Disk Formation in Massive Galaxies at Cosmic Dawn

Introduction

The James Webb Space Telescope (JWST) has made it possible to measure the sizes and internal structure of galaxies during the first few hundred million years of cosmic history. Galaxy sizes encode how stellar mass is distributed within the gravitational potential (for reviews, see Conselice 2014; Naab & Ostriker 2017), and thus constrain the angular momentum of the accreted gas (Fall & Efstathiou 1980; Mo et al. 1998; Dutton & van den Bosch 2012; Somerville et al. 2018), the spatial distribution of star formation (Nelson et al. 2016), and the strength of feedback (Sales et al. 2010; El-Badry et al. 2016).

At fixed luminosity or stellar mass, galaxies become progressively smaller towards higher redshift, with effective radii scaling roughly as $(1+z)^{-1}$ over $z \sim0$–10 (Shibuya et al. 2015; Yang et al. 2025). JWST has extended these measurements to $z > 9$, where galaxies have UV half-light radii of a few hundred parsecs (Ono et al. 2023, 2025; Morishita et al. 2024, hereafter M24). Over $5 < z < 14$, however, M24 find a much slower average size evolution, $R_{\rm e} \propto(1+z)^{-0.4\pm0.2}$, suggesting that the decline in size flattens at cosmic dawn. This average conceals a compact population: $\sim13$ per cent of their sample is only marginally resolved in NIRCam imaging, with $R_{\rm e} \lesssim100$ pc and star-formation rate (SFR) surface densities of $\Sigma_{\rm SFR} > 10\ M_\odot\,{\rm yr}^{-1}\,{\rm kpc}^{-2}$, a regime occupied by fewer than 0.01 per cent of local star-forming galaxies.

Subsequent spectroscopic studies have examined this compact population more systematically and suggest that the size distribution of luminous galaxies at $z > 10$ may be bimodal (Harikane et al. 2025; Roberts-Borsani et al. 2026). The most compact of these sources preferentially show strong high-ionization UV lines, more negative UV continuum slope $\beta$, and high SFR surface densities (Harikane et al. 2025; Roberts-Borsani et al. 2026). Rather than treating them as a distinct class, Roberts-Borsani et al. 2026 interpret these compact strong-line sources as a short-lived phase near the peak of a central starburst along a common evolutionary pathway, during which galaxies temporarily occupy the compact end of the size distribution.

The size–mass relation offers a complementary constraint on structural evolution. M24 find a positive UV size–mass relation at $5 < z < 14$, with the compact systems at $M_\star< 10^9\,M_\odot$ forming a low-mass tail that lies $\gtrsim1.5\sigma$ below the main relation. These compact sources are more prevalent in their lensing-cluster field, where magnification aids the detection of small, low-mass galaxies. Lensed faint galaxies can be compact down to tens of parsecs, comparable to individual star-cluster complexes (Bouwens et al. 2017; Adamo et al. 2024; Mowla et al. 2024), and their size–luminosity relation is much steeper than the normal one (Bouwens et al. 2017). In the rest-frame optical, JWST samples show a positive size–mass relation at $2 < z < 8$ with a non-evolving slope of $\simeq0.2$–$0.25$ (Miller et al. 2025; Yang et al. 2025). In addition, a population of blue, barely resolved sources accounts for roughly 10 per cent of the sample at $z \gtrsim5$–$6$ (Miller et al. 2025; Yang et al. 2025). At the high-mass end, Chen et al. (2026b) find that the optical size–mass relation flattens above a nearly redshift-independent stellar mass of $\sim10^{10}\,M_\odot$ over $0.5 < z < 6$, and the bending becomes progressively stronger at higher redshifts.

Cosmological simulations diverge in their predicted size–mass relation at high redshift, as well as in their treatment of compact objects at the low-mass end. A positive stellar half-mass radius–mass relation has been found in FIRE-2 at $z = 5$–$10$ ([85]), in the SERRA galaxies at $z \simeq7.7$ (Pallottini et al. 2022), and in THESAN-ZOOM at $z = 3$–$13$ (McClymont et al. 2025). Of these, the corresponding UV half-light size–mass relation in FIRE-2 and THESAN-ZOOM are closer to the observed one for typical galaxies, while SERRA agrees better with lensed galaxies at $6 \leq z \leq8$. Moreover, many other simulations have difficulty reproducing a positive size–mass relation at high redshift, let alone the compact galaxies at the low-mass end. A negative stellar half-mass radius–mass relation above $M_\star\simeq5 \times10^8\,M_\odot$ has been reported in IllustrisTNG (Genel et al. 2018), BLUETIDES (Marshall et al. 2022), FLARES (Roper et al. 2022), THESAN-1 (Shen et al. 2024), and FirstLight (Cataldi et al. 2026). These studies commonly attribute the negative relation to wet compaction at high redshift, although a positive UV half-light radius–mass relation can be recovered once dust effects are included, as in THESAN-1 (Shen et al. 2024) and FirstLight (Ceverino et al. 2026).

Wet compaction connects gas inflow, central star formation, and structural evolution, and is widely used to explain the negative trend in the size–mass relation at high redshift. Dekel & Burkert (2014) first proposed that gas-rich disks can contract into compact star-forming “blue nuggets” when the inflow time is shorter than the star-formation time, allowing gas to reach the centre before being converted into stars. Cosmological simulations by Zolotov et al. (2015) trace this dissipative buildup and the subsequent central gas depletion and inside-out quenching, with lower-mass systems experiencing repeated compaction episodes before sustained quenching. Tacchella et al. (2016) further show how cycles of gas compaction, depletion, and replenishment produce excursions around the star-forming main sequence. Dekel et al. (2019) link the preferred mass scale of major compaction to the combined effects of supernova feedback in low-mass haloes and virial shock heating in more massive haloes. Lapiner et al. (2023) identify major compaction events near $M_\star\sim10^{10}\,M_\odot$, accompanied by a transition from dark-matter to baryon dominance in the centre. These compaction events can be triggered by wet mergers, counter-rotating streams, violent disk instability, and gas recycling (Zolotov et al. 2015; Dekel et al. 2019; Lapiner et al. 2023).

Compaction can promote disk formation and inside-out quenching. In cosmological simulations, gas depletion after the central starburst suppresses star formation in the compact core, while an extended star-forming disk can grow around it (Zolotov et al. 2015; Tacchella et al. 2016; Lapiner et al. 2023). In TNG50, the compact centres that compaction builds in massive galaxies are themselves flattened and rotation supported, resembling nuclear disks and bars rather than classical bulges, and extended disks assemble around them before the centres become dispersion dominated at $z \lesssim1$ (Luo et al. 2026). Controlled simulations and analytical models suggest that a concentrated central potential can help gas settle into coherent rotation (Hopkins et al. 2023; Semenov 2026). These results suggest that size evolution, disk formation, and inside-out quenching may be coupled through central stellar buildup and the subsequent redistribution of star formation.

Observations also provide evidence for a connection between compact central structures and rotation at early epochs. Baker et al. (2025) find a dense core surrounded by an extended star-forming component at $z = 7.43$, with a radial increase in specific SFR indicating inside-out growth. Resolved gas kinematic maps establish a dynamically cold disk in REBELS-25 at $z = 7.31$ (Rowland et al. 2024, 2026), while the Cosmic Grapes galaxy at $z = 6.072$ shows that numerous star-forming clumps can coexist with a rotating gas disk (Fujimoto et al. 2025). At $z = 10.6$, recent spectroscopy of GN-z11 supports a substantial contribution from young massive stars to its compact emission (Chen et al. 2026a; Nakane & Ouchi 2026), and its kinematics show tentative rotational support with $v/\sigma\simeq3$ ([151]). More indirectly, in a photometric sample of 169 candidates at $z \sim10$–$16$, Ono et al. (2025) find a log-normal size distribution, a nearly uniform distribution of projected axis ratios over $q = 0.2$–$0.8$, approximately exponential profiles with Sérsic indices close to unity, and a characteristic size-to-halo-radius ratio of $R_{\rm e}/R_{\rm v} \simeq0.015$, similar to the nearly constant ratio of $\simeq0.01$–$0.035$ found over $z \sim0$–$8$ (Shibuya et al. 2015; Yang et al. 2025). These properties are consistent with the classical disk-formation picture, in which disk sizes are set by the angular momentum that galaxies inherit from their haloes and therefore scale with the halo virial radius (Fall & Efstathiou 1980; Mo et al. 1998), although morphology alone does not establish rotational support. Together, these results motivate a study of central buildup and disk formation beyond $z = 9$, where it remains an open question whether this connection is already in place.

In this paper, we introduce MAGE (Massive Assembly of Galaxies in the Early universe), a new suite of 11 cosmological zoom-in simulations run with the RAMSES adaptive mesh refinement (AMR) cosmological hydrodynamics code, to examine these connections over $9 \lesssim z \lesssim15$. With $\sim10$ pc resolution and the star-formation and feedback model of Andalman et al. (2025), the simulations allow us to resolve compact galaxies more reliably and to study their connection to disk formation at $z > 9$.

This paper is organized as follows. § 2 describes the simulations and analysis methods. § 3 presents the size evolution of representative galaxies and the statistical size–mass and size–magnitude relations. § 4 examines in-situ mass fractions and the spatial distributions of stellar ages and star formation. § 5 connects the emergence of rotational support to central mass buildup and gas flows. § 6 discusses wet compaction, the conditional disk-formation diagnostic, feedback dependence, and comparisons with simulations and observations. We summarize the results in § 7.

Methods

Numerical simulations

The MAGE suite comprises 11 cosmological zoom-in simulations of massive haloes evolved to $z \simeq9$ with RAMSES. The suite extends the fiducial simulation of Andalman et al. (2025) using the same star-formation and feedback prescriptions. We summarize the numerical methods here and refer to that work for implementation details.

We assume a Planck cosmology ($\Omega_{\mathrm{m}} = 0.308$, $\Omega_{\Lambda} = 0.692$, $\Omega_{\mathrm{b}} = 0.0484$, $H_{0} = 67.81\ \mathrm{km\ s^{-1}\ Mpc^{-1}}$; Planck Collaboration et al. [116]), except for MAGE07, the fiducial run of Andalman et al. (2025), which adopts slightly different cosmological parameters ($\Omega_{\mathrm{m}} = 0.276$, $\Omega_{\Lambda} = 0.724$, $\Omega_{\mathrm{b}} = 0.049$, $H_{0} = 70\ \mathrm{km\ s^{-1}\ Mpc^{-1}}$). The target haloes are selected from large parent boxes of size $100 h_{0}^{-1}\ \mathrm{cMpc}$ and $250 h_{0}^{-1}\ \mathrm{cMpc}$ ($h_{0} \equiv H_{0}/(100\ \mathrm{km\ s^{-1}\ Mpc^{-1}})$), chosen to host haloes with $M_{\mathrm{v}} > 10^{10}\ M_{\odot}$ by $z = 9$. To construct zoom-in initial conditions, we trace the dark matter particles within three times the virial radius of each target halo back to their positions at $z \simeq124$, the starting redshift of the parent boxes, and define the convex hull enclosing these particles as the zoom region using MUSIC (Hahn & Abel (2011)). The zoom-in initial conditions are generated at $z = 50$, and the zoom-in simulations are started from that epoch. The zoom-in simulations are run with the AMR code RAMSES (Teyssier (2002)), which solves the Euler equations with a second-order Godunov scheme, the HLLC Riemann solver, and MUSCL–Hancock reconstruction. Dark matter and star particles are evolved in a Lagrangian framework, while gas is modelled on an Eulerian grid that refines dynamically according to mass thresholds and local Jeans criteria. This approach yields an effective physical resolution of $\sim10\ \mathrm{pc}$ at the highest refinement level throughout the simulations. The main properties of each zoom-in halo are summarized in Table 1.

Name	$L_{\mathrm{box}}$	$m_{\mathrm{dm,min}}$	$m_{\mathrm{b,min}}$	$\Delta x_{\mathrm{min}}$	$M_{\mathrm{v}}$	$M_{\star}$	$M_{\mathrm{g}}$	$R_{\mathrm{v}}$	$R_{\mathrm{e}}$	$R_{1/2}$	$\mathrm{SFR}_{50}$	$(v/\sigma)_{\mathrm{SFgas}}$	$(v/\sigma)_{\star,10}$	$(v/\sigma)_{\star}$
MAGE01	250	4.72	8.80	8.80	6.36	31.56	7.80	26.70	0.16	0.22	209.47	1.35	3.52	0.80
MAGE02	250	4.72	8.80	8.80	3.45	14.03	6.73	21.76	0.26	0.53	119.01	1.17	1.50	0.78
MAGE03	250	4.72	8.80	8.80	3.14	13.12	5.74	21.14	0.23	0.33	104.69	3.50	3.20	1.00
MAGE04	100	2.42	4.50	7.03	2.23	10.93	3.46	18.84	0.21	0.27	62.10	1.70	2.50	1.10
MAGE05	250	4.72	8.80	8.80	2.10	3.46	2.59	18.45	0.19	0.31	31.38	1.08	1.37	0.56
MAGE06	250	4.72	8.80	8.80	1.33	2.92	2.48	15.89	0.22	0.34	18.16	2.00	2.40	1.59
MAGE07	100	2.04	4.36	6.81	0.95	3.67	3.13	14.17	0.18	0.24	45.90	1.26	1.69	1.16
MAGE08	250	4.72	8.80	8.80	0.44	0.99	0.61	11.08	0.04	0.10	2.19	2.06	1.13	0.46
MAGE09	250	4.72	8.80	8.80	0.36	0.38	0.39	10.28	0.11	0.21	2.83	0.87	1.60	0.64
MAGE10	250	4.72	8.80	8.80	0.27	0.50	0.35	9.47	0.04	0.06	2.31	0.47	0.69	0.83
MAGE11	250	4.72	8.80	8.80	0.27	0.23	0.07	9.42	0.02	0.03	0.89	0.57	1.25	0.42

To capture unresolved turbulent motions, we employ a Large Eddy Simulation (LES) framework that evolves the turbulent kinetic energy (TKE) using a sub-grid model based on Kretschmer & Teyssier (2020). The sub-grid model decomposes density into a volume-averaged component and fluctuations, while velocity and temperature are decomposed using mass-weighted (Favre) averages (Schmidt & Federrath (2011)). The sub-grid turbulence model tracks the evolution of the TKE with creation and dissipation source terms determined from the local velocity gradients and a dissipation time-scale set by the grid scale and turbulent velocity. Turbulent pressure contributes to supporting gas against collapse and is particularly relevant in regulating star-formation efficiency (SFE) in massive, rapidly accreting systems at high redshift.

Star formation proceeds according to a multi-freefall model (Federrath & Klessen (2012)) as implemented in RAMSES by Kretschmer & Teyssier (2020). The local star-formation rate density follows the Schmidt law, $\dot{\rho}_{\star} = \epsilon_{\mathrm{ff}}\rho_{\mathrm{g}}/t_{\mathrm{ff}}$, with the star-formation efficiency per free-fall time $\epsilon_{\mathrm{ff}}$ depending on the local gas properties and sub-grid turbulence. In contrast to conventional fixed-efficiency models, this framework allows $\epsilon_{\mathrm{ff}}$ to vary dynamically in response to the turbulent state of the gas. The fiducial model assumes fully compressive turbulence forcing, which leads to elevated $\epsilon_{\mathrm{ff}}$ in dense regions.

Supernova (SN) feedback is implemented via a hybrid thermal and kinetic injection model. Each SN injects thermal energy locally and radial momentum through the neighbouring cell faces, with the momentum calibrated using cooling-radius scaling relations from Martizzi et al. (2015). The feedback is purely local and does not resolve the Sedov–Taylor expansion in detail, but approximates the net momentum transfer over unresolved stages of remnant evolution. The number of SNe per star particle is set by two parameters: the fraction of each star particle’s mass in SN progenitors, 0.2, and the mean progenitor mass, $10\,M_{\odot}$. Each solar mass of stars formed therefore produces $0.2/10 = 0.02$ SNe. Metal enrichment is included with a yield of $y = 0.1$, so each SN returns 10 per cent of the progenitor mass, $1\,M_{\odot}$, as metals.

In addition to SN feedback, an early feedback model accounts for photoionization heating from massive stars. Gas in cells containing young star particles (age $< 20\ \mathrm{Myr}$) is heated to $T_{\mathrm{phot}} = 2 \times10^{4}\ \mathrm{K}$ to mimic the thermal pressure of ionized H II regions. This heating is implemented by setting the thermal energy corresponding to a sound speed of $c_{\mathrm{s}} \simeq22\ \mathrm{km\ s^{-1}}$, assuming an ideal monatomic gas. While this model neglects radiation pressure and stellar winds (see § 6.8), it captures the rapid thermal response of star-forming gas to early stellar feedback.

Gas cooling is computed using equilibrium chemistry for hydrogen and helium (Katz et al. (1996)), with a metallicity-dependent cooling curve from Sutherland & Dopita (1993). A uniform, time-evolving UV background (Haardt & Madau (1996)) is switched on at $z = 10$, with self-shielding that suppresses it in gas denser than $n_{\mathrm{H}} \sim0.01\ \mathrm{cm^{-3}}$. It is therefore active only during the final $\simeq70\ \mathrm{Myr}$ of our analysis ($9 \lesssim z < 10$). We find no discontinuity at $z = 10$ in the gas content of the galaxies, and the background is expected to have little effect on the star-forming gas over this short interval. A temperature floor is imposed at the redshift-dependent cosmic microwave background (CMB) temperature as an approximate treatment of heating by CMB photons. The initial gas metallicity is set to $Z = 10^{-3}Z_{\odot}$, so the simulations model enriched, Population II star formation without explicitly following Population III enrichment. Limitations of the sub-grid prescriptions are discussed by Andalman et al. (2025).

Halo finders, merger trees, and galaxy properties

Halo catalogues and merger histories are constructed with the six-dimensional phase-space halo finder ROCKSTAR (Behroozi et al. 2013a), coupled with the gravitationally consistent merger tree algorithm CONSISTENT TREES (Behroozi et al. 2013b). The target halo is traced along the main branch of its merger tree, which extends to $z \sim20$–$24$.

Each simulation stores 27–49 outputs (426 in total) between the initial conditions at $z = 50$ and the final output at $z \simeq9$. Only a few outputs are written above $z \simeq16$, whereas at $z \lesssim16$ they are spaced by a median of $\simeq10$ Myr (16th–84th percentiles 6–13 Myr). We analyse the target galaxies only from $z \sim15$–16 onwards: at earlier times the progenitors have masses of just $\sim10^{6}$–$10^{7}\,M_{\odot}$ and the outputs are sparse, so halo and galaxy properties are not robust. This yields 20–39 analysed snapshots per target galaxy and 318 in total.

The halo catalogues are constructed with ROCKSTAR applied to the dark matter particles alone, so the virial mass and radius it returns describe only the dark matter component of the hydrodynamical simulations. We therefore recompute both quantities from the total matter distribution. Around the ROCKSTAR centre of each distinct (host) halo, we sum the enclosed mass of dark matter particles, star particles and gas cells as a function of radius. We define the virial radius $R_{\rm v}$ as the radius at which the mean enclosed total density first falls below $\Delta_{\rm vir}(z)\,\rho_{\rm crit}(z)$, where $\Delta_{\rm vir}$ follows Bryan & Norman (1998), the same virial definition adopted by ROCKSTAR. Because $\Omega_{\rm m}(z) \simeq1$ at the redshifts considered here, this threshold corresponds to $\simeq178$ times the mean matter density. The virial mass is then $M_{\rm v} = (4\pi/3)\,\Delta_{\rm vir}\rho_{\rm crit}R_{\rm v}^{3}$, i.e. all dark matter, stars and gas within $R_{\rm v}$, so that the virial radius and mass are defined self-consistently.

The host and subhalo relationships are taken from CONSISTENT TREES. Merger mass ratios are defined as the ratio of the infalling halo mass to the host halo mass, evaluated at the final snapshot before the infalling halo becomes a subhalo. This minimizes the impact of subsequent tidal stripping.

For the statistical sample, the galaxy centre is determined with an iterative shrinking-sphere method (Power et al. 2003). Starting from a sphere of radius $0.3\,R_{\rm v}$ centred on the halo centre returned by ROCKSTAR, we compute the mass-weighted centre of the enclosed star particles, recentre the sphere on this point, and shrink its radius by a factor of 0.7 at each iteration. We repeat this procedure down to a radius of $0.02\,R_{\rm v}$, terminating early if fewer than 20 star particles remain enclosed, and adopt the final mass-weighted stellar centre as the galaxy centre. This procedure reduces the influence of satellites and offsets between the stellar and halo centres. Visual inspection nevertheless reveals occasional miscentring, particularly during major mergers when multiple dense stellar components lie within $0.3\,R_{\rm v}$. For the 11 target central galaxies, we visually inspect the centres and manually fix selected snapshots to follow the main stellar component consistently. These corrections affect only $\sim5$ per cent of the snapshots, so we expect the statistical results to be largely unaffected.

We define the galaxy boundary as $R_{\rm gal} \equiv0.15\,R_{\rm v}$. Unless a different aperture is specified, stellar mass, gas mass, SFR, and the 3D stellar half-mass radius $R_{1/2}$ are computed within this boundary at each snapshot. The enclosed mass, $M(<r)$, is the cumulative mass within a sphere of radius $r$. We define the circular speed as $v_{\rm c}(r) = \sqrt{GM(<r)/r}$ and the halo virial speed as $v_{\rm v} = \sqrt{GM_{\rm v}/R_{\rm v}}$. Component circular speeds use the corresponding enclosed mass. We also tested a larger galaxy boundary of $R_{\rm gal} = 0.3\,R_{\rm v}$. The derived stellar properties are largely insensitive to this choice, whereas the gas properties exhibit a markedly stronger dependence on the adopted boundary (see also Stevens et al. 2014).

Gas inflow and outflow rates across $R_{1/2}$, $R_{\rm gal}$, and $R_{\rm v}$ are estimated by multiplying the volume average of $\rho v_{\rm r}$ in a spherical shell of thickness $\Delta R = 0.1\,R$ centred on each radius $R$ by the corresponding spherical area, $4\pi R^{2}$. Inward- and outward-moving gas are treated separately, and both rates are plotted as positive magnitudes. Radial velocities are measured relative to the mass-weighted stellar bulk velocity within $R_{1/2}$. SFRs are estimated by summing the formed masses of stars younger than the adopted 5 or 50 Myr window within the present aperture and dividing by that window. If $m_i$ is its current mass and $t_i$ its age, we adopt $m_{\rm form,i} = m_i/f_{\rm SN}(t_i)$, where

is the mean fraction of its initial mass that a particle retains under our SN feedback model, which returns 20 per cent of it to the gas at a constant rate between 3 and 20 Myr. This is the same estimate RAMSES uses to set each particle’s SN rate. Because individual SNe are drawn stochastically, a particle’s true initial mass differs from $m_{\rm form,i}$ by $\lesssim2$ per cent, which is negligible once summed over a galaxy. The 5 Myr average emphasizes recent bursts and overlaps the onset of SN feedback, whereas the 50 Myr average smooths the burst–lull cycle. The two intervals also correspond to the time-scales probed by nebular emission lines and by the rest-frame UV continuum, respectively (Kennicutt & Evans 2012; Flores Velázquez et al. [54]), so their ratio is the simulation analogue of the H$\alpha$-to-UV burstiness indicator used at high redshift (Faisst & Morishita 2024; Dome et al. 2024).

We classify stars by their birth positions relative to the main progenitor: a star is in situ if it formed within the galaxy aperture ($0.15R_{\rm v}$) at its birth time, and ex situ otherwise. This radial criterion misclassifies stars formed in companion galaxies that have already entered the aperture during mergers. To correct for this, we identify spatially distinct, persistent star-forming companions within the aperture by applying the density-based clustering of star particles with noise (DBSCAN) algorithm (Ester et al. 1996) to the birth positions of stars formed in short ($\simeq2$ Myr) time windows, short enough that a companion does not move appreciably within a window. The main galaxy is the cluster dominating the central few hundred parsecs, and stars born in other sufficiently massive clusters that persist over consecutive windows are reclassified as ex situ. Small star-forming clumps within the disk are thus not counted as companions. This step can only move stars from in situ to ex situ. Its effect is largest during mergers, but at $z = 9$ it reduces the in-situ fraction by a median of 0.04 and changes the half-mass radius of the in-situ population by at most 10 per cent.

For each in-situ star we record its birth radius, $R_{\rm birth}$, defined as the distance between its birth position and the centre of the main progenitor at its formation time, $t_{\rm birth}$, together with the properties of the main progenitor at that epoch: its stellar mass within $R_{\rm gal}$, $M_{\star}(t_{\mathrm{birth}})$, its stellar half-mass radius, $R_{1/2}(t_{\mathrm{birth}})$, and its virial radius, $R_{\mathrm{v}}(t_{\mathrm{birth}})$. Birth positions are stored by the simulation at the time of formation. Because the outputs are spaced by $\simeq10\,\mathrm{Myr}$, the progenitor centre and properties at $t_{\mathrm{birth}}$ are linearly interpolated in time between the two bracketing outputs.

We employ all galaxies located within the zoom-in region for our statistical analysis, provided that their host haloes are distinct in the CONSISTENT TREES merger trees. Specifically, in addition to the principal target central galaxy of each zoom-in simulation, we include all galaxies that satisfy the following criteria:

1. $\log_{10}(M_{\star}/M_{\odot}) \gtrsim7.5$;

2. $9 \lesssim z \lesssim15$;

3. central galaxies of distinct haloes;

4. a fraction of high-resolution dark-matter particles within $R_{\mathrm{v}}$ such that $N_{\mathrm{dm,highres}}/N_{\mathrm{dm,tot}} \geq0.9$;

5. a fraction of the total stellar mass beyond half the galactic radius satisfying $f(R > R_{\mathrm{gal}}/2) < 0.25$.

A central galaxy in criterion (iii) is the galaxy located within $R_{\mathrm{gal}} = 0.15\,R_{\mathrm{v}}$ of the centre of a distinct host halo in the merger trees. Galaxies hosted by subhaloes are treated as satellites and excluded. The high-resolution particles in criterion (iv) have $m_{\mathrm{dm}} < 5 \times10^{4}\,M_{\odot}$. The stellar-mass threshold in criterion (i) selects systems with more than 3000 star particles within $R_{\mathrm{gal}}$. Below this threshold, measured sizes approach the numerical grid scale and the gravitational softening length of star particles, which is linked to the local AMR grid in RAMSES. The threshold of 0.25 in criterion (v) is motivated by the conventional 1:4 stellar-mass ratio used to define major mergers ((Puskás et al. 2025)). We adopt this value as an empirical cut on the outer stellar-mass fraction after testing different thresholds and inspecting their effectiveness in removing systems with multiple spatially separated stellar components. Our relatively large aperture, $R_{\mathrm{gal}} = 0.15R_{\mathrm{v}}$, places $R_{\mathrm{gal}}/2$ well beyond the typical stellar half-mass radius, so the cut removes multi-component systems without penalizing extended single galaxies. Excluding criterion (v) only increases the upper scatter of $R_{1/2}$ and does not significantly affect the median in any stellar-mass or redshift bin.

The selection yields 2415 galaxy–snapshot measurements over the analysed redshift range. These are repeated measurements of evolving systems, rather than 2415 independent galaxies. Their distributions in the halo-mass–redshift and stellar-mass–redshift planes are shown in Figure 1. The non-target galaxies primarily populate the low-mass regime, whereas the high-mass end remains dominated by the target central galaxies, marked by red stars at $z \simeq9$. For statistics presented as a function of stellar mass, we divide the sample into three redshift intervals over $9 \lesssim z \lesssim15$ and bin the stellar mass with a width of $\Delta\log_{10}(M_{\star}/M_{\odot}) = 0.5$. For statistics presented as a function of mass, we divide the sample into three stellar-mass intervals over $7.5 \leq\log_{10}(M_{\star}/M_{\odot}) \leq10.5$ and bin redshift with a width of $\Delta z = 0.5$. We show bins containing at least seven galaxy–snapshot measurements. Unless otherwise specified, curves show medians and error bars span the 16th–84th percentiles of these measurements.

Image not extracted

For reference, the middle panel of Figure 1 places the sample relative to theoretical thresholds for efficient star formation at cosmic dawn. Feedback-free starbursts (FFBs) are expected in haloes above $M_{\mathrm{v}} \simeq10^{10.8}[(1+z)/10]^{-6.2}\,M_{\odot}$, where gas densities exceed $\sim3 \times10^{3}\,\mathrm{cm}^{-3}$ and free-fall times become shorter than the $\sim1\,\mathrm{Myr}$ delay before stellar winds and SNe act ((Dekel et al. 2023; Li et al. 2024)). This threshold nearly coincides with the halo mass above which central dark-matter accelerations exceed the critical value at which stellar feedback can no longer overcome gravity, $M_{\mathrm{v}} \simeq6.5 \times10^{10}[(1+z)/10]^{-6}\,M_{\odot}$ ((Boylan-Kolchin 2025)). In a

closely related formulation, (Somerville et al. 2025) let the cloud-scale star-formation efficiency rise smoothly with gas surface density, approaching unity above $\Sigma_{\mathrm{crit}} \simeq2 \times10^{3}\,M_{\odot}\,\mathrm{pc}^{-2}$, where gravity overcomes the momentum injected by stellar winds and radiation pressure (see also Grudić et al. 2018; Menon et al. 2024). Implemented in a semi-analytic model, this prescription raises the predicted abun­ dance of UV-luminous galaxies at $z \simeq12$–14 by one to two orders of magnitude and reproduces the observed UV luminosity functions. At $z = 9$, the haloes of MAGE01–MAGE07 ($M_{\rm v} \simeq1$–$6 \times10^{11}\,M_\odot$, see Table 1) lie above both halo-mass thresholds, whereas those of MAGE08–MAGE11 ($M_{\rm v} \simeq3$–$4 \times10^{10}\,M_\odot$) lie below them but above the halo mass at which the virial velocity reaches the characteristic SN-feedback scale $v_{\rm v} \simeq100\,\mathrm{km\,s^{-1}}$ (Dekel & Silk 1986, see § 5). The middle panel also shows contours of constant comoving number density, $n(>M_{\rm v})$, from the virial halo mass function of Despali et al. (2016). The bulk of the sample lies at halo masses with $n(>M_{\rm v}) \gtrsim10^{-3}\,\mathrm{Mpc}^{-3}$, whereas the target haloes at $z \simeq9$ are rare objects with $n(>M_{\rm v}) \simeq10^{-7.6}$–$10^{-3.4}\,\mathrm{Mpc}^{-3}$.

The bottom panel of Figure 1 marks the characteristic stellar mass of major wet-compaction events in cosmological simulations at cosmic noon, $M_\star\simeq10^{10}\,M_\odot$ (Lapiner et al. 2023), a scale we return to in § 6.1. Only the four most massive targets (MAGE01–MAGE04) reach this mass by $z \simeq9$, and galaxy–snapshot measurements above it are confined to $z \lesssim10.5$. The bulk of the sample lies one to two decades below it, so the sample crosses the compaction scale only in its high-mass tail.

Kinematics of star-forming gas and stars

We characterize rotational support with the rotation-to-dispersion ratio, $v/\sigma$ (often denoted $v_{\rm rot}/\sigma$ or $V/\sigma$ in the literature), measured within $R_{\rm gal} = 0.15R_{\rm v}$ for three tracers: the star-forming gas, stars younger than 10 Myr, and the full stellar population. All three are measured with the same procedure, which differs only in the weight $w$ assigned to each resolution element (the SFR for gas cells and the mass for star particles) and, for the gas, an additional pressure term. Each tracer defines its own frame from its $w$-weighted bulk velocity and total angular momentum. After subtracting the bulk velocity, we align the angular-momentum vector with the $z$-axis and calculate the cylindrical velocity components $v_R$, $v_\phi$, and $v_z$. The rotation speed and effective one-dimensional velocity dispersion are

where $\langle\cdot\rangle_w$ denotes the $w$-weighted mean, $\sigma_k^2 = \langle(v_k-\langle v_k\rangle_w)^2\rangle_w$ is the variance of the resolved velocity component $v_k$, and $\sigma_P$ accounts for pressure support that is not carried by resolved motions. For stars, $\sigma_P = 0$.

For the star-forming gas, SFR weighting selects the gas that is forming stars at that instant, the component traced observationally by nebular lines. Because this gas is supported by thermal and unresolved turbulent pressure in addition to its resolved bulk motions, we set $\sigma_P^2 = \langle c_s^2+\sigma_{\rm turb}^2\rangle_w$, where $c_s$ is the thermal sound speed and $\sigma_{\rm turb}$ the one-dimensional sub-grid turbulent velocity of the LES model (§ 2). We denote the resulting quantities $v_{\rm rot,SFGas}$ and $\sigma_{\rm SFGas}$, and their ratio $(v/\sigma)_{\rm SFGas}$.

The two stellar tracers bracket the kinematic memory of the stellar population. Young stars trace the recently formed stellar component and initially inherit the resolved velocities of their parent gas, therefore they are the collisionless counterpart of the star-forming gas, without the pressure term. The full population spans many orbits and several burst–lull cycles, and tests whether rotational support survives phase mixing and feedback. For our simulated tracers, $v/\sigma$ is defined by Eq. (2), whereas quoted observational values retain their authors’ tracers and dispersion definitions.

Rest-frame UV surface brightness

For each selected central galaxy, we construct intrinsic (dust-free, noiseless) JWST/NIRCam F200W images using the Flexible Stellar Population Synthesis code (FSPS; Conroy et al. 2009; Conroy & Gunn 2010). We adopt simple stellar populations (SSPs) with a Chabrier initial mass function (Chabrier 2003). MIST isochrones, and the MILES stellar spectral library. Nebular emission and internal dust attenuation are omitted. F200W samples rest-frame UV wavelengths of approximately 2000 Å at $z = 9$ and 1250 Å at $z = 15$. Photometry is tabulated over stellar age and metallicity, with each particle’s metallicity mapped to the nearest available FSPS index. The SSP flux density, normalized per unit formed stellar mass, is scaled by the particle’s birth mass per unit formed stellar mass, $m_{\mathrm{form},i}/M_\odot$. Thus the particle flux weight is $f_{\nu,i}=f_{\nu,\mathrm{SSP}}(t_i,Z_i,z)\,m_{\mathrm{form},i}/M_\odot$, with the SSP flux obtained from its AB magnitude (Oke & Gunn 1983).

We select stars within the three-dimensional $R_{\rm gal}$ sphere and sum their flux weights onto a $256 \times256$ grid spanning $-R_{\rm gal}$ to $+R_{\rm gal}$ in each projected coordinate. The measurement field therefore has width $2R_{\rm gal}$ and pixel size $2R_{\rm gal}/256$, which ranges from a few parsecs for the least massive haloes to $\simeq30$ pc for the most massive target halo at $z = 9$, comparable to the $\sim10$ pc resolution of the simulations. It is several to tens of times finer than the NIRCam short-wavelength pixel scale of $0.031''$ ($\simeq130$ pc at $z = 10$) and far below the F200W point spread function (PSF) full width at half-maximum (FWHM) of $\simeq0.066''$ ($\simeq280$ pc). These images are therefore intrinsic surface-brightness maps, intended for illustration and for comparing intrinsic sizes among simulations, rather than mock observations. Reproducing NIRCam measurements would require dust attenuation, PSF convolution, noise, observational pixelization, and the source-detection and profile-fitting procedures of the observations. We will present mock observations in more detail, including the effects of dust, in a follow-up study (Nakazato et al., in preparation).

UV half-light and half-mass radii

We measure intrinsic UV sizes from projected F200W images, which include neither noise nor PSF convolution. Sources are identified with PHOTUTILS (Bradley et al. 2024) as regions of at least five connected pixels above a threshold of one standard deviation above the median of the image. Because the images are noiseless, this threshold acts as a surface-brightness contour that delineates individual galaxies and clumps rather than as a signal-to-noise cut. For the target-galaxy histories we select the source whose centroid is closest to the image centre, and for the statistical catalogue the source with the largest area. Absolute UV magnitudes are calculated from the FSPS rest-frame spectrum near 1500 Å, summing the contributions of all particles assigned to the same selected source.

We define the UV half-light radius, $R_{\rm e}$, and the projected stellar half-mass radius, $R_{1/2,\mathrm{2D}}$, using the same non-parametric estimator, the circular half-light radius (rhalf_circ) from STATMORPH (Rodriguez-Gomez et al. 2019), applied to the F200W image and to a projected stellar-mass image constructed on the same pixel grid, respectively. In each case it is the radius of the circle, centred on the asymmetry-minimizing centre of the image, that encloses half of the light (or mass) within the outer radius rmax_circ. Both measurements use the source and segmentation map defined from the F200W image: emission from other detected sources is excluded, while faint emission below the detection threshold is included. The two radii therefore differ only in their weighting by UV light or stellar mass.

Our primary intrinsic size is the three-dimensional stellar half-mass radius, $R_{1/2}$: the radius of the sphere, centred on the galaxy centre of § 2.2, that encloses half of the stellar mass within $R_{\rm gal}$. Unlike $R_{\rm e}$ and $R_{1/2,2{\rm D}}$, it depends on neither the viewing direction nor the source segmentation. The three radii thus form a sequence: $R_{1/2}$ measures the intrinsic mass distribution, $R_{1/2,2{\rm D}}$ the projected mass distribution of the segmented source, and $R_{\rm e}$ the projected UV light of the same source.

Size Evolution Across Cosmic Time

We first compare the size histories of the smooth-accretion system MAGE07 and the major-merger system MAGE04 (§ 3.1). We then examine the size–mass and size–redshift relations of the full statistical sample (§ 3.2), followed by the UV size–magnitude relation and its observational comparison (§ 3.3).

Size evolution of two example galaxies

MAGE07 and MAGE04 illustrate contrasting assembly histories at cosmic dawn. The former grows mainly through smooth accretion, whereas the latter has undergone recent major mergers. Despite these differences, they reach similar rest-frame UV morphologies and sizes by $z \simeq9$.

Figure 2 traces the evolution of MAGE07 through its rest-frame UV morphology, mass buildup, mass flow rates, and characteristic size. The top panels present intrinsic F200W-band surface-brightness maps at $z = 12.37$, $z = 10.55$, and $z = 9$. At the earliest time, the emitting region spans only a few image pixels, with most of the UV light concentrated near the centre. As the system evolves, the emitting region becomes progressively more extended and develops multiple star-forming clumps within the galaxy. By $z = 9$, MAGE07 has transformed from a compact progenitor into an extended galaxy.

Image not extracted

The lower panels of Figure 2 show that this morphological growth is accompanied by a relatively smooth buildup of mass. Halo mergers onto the main branch with mass ratios greater than 0.01 are indicated by lollipops at the top of each panel. No major halo merger occurs for this system at $z < 14$, although a major galaxy merger at $z \simeq13.5$ leaves a visible imprint, as discussed below. The halo mass therefore grows steadily, while the gas and stellar masses within $R_{\rm gal}$ increase overall with time. By $z = 9$, the stellar mass reaches $M_\star\simeq0.22\,f_{\rm b}M_{\rm v}$, where $f_{\rm b}$ is the cosmological baryon fraction, slightly above the uppermost grey line, which marks a global star-formation efficiency of 0.2. The mild, transient dip in stellar mass at $t \simeq385$ Myr is caused by orbiting satellites temporarily moving outside $R_{\rm gal}$ at this snapshot. The rapid consumption of gas at $t \simeq325$ Myr is driven by a major galaxy merger, an event accompanied by a transient peak in $R_{1/2}$. Despite this depletion, the system remains gas-rich throughout most of its evolution, maintaining a gas mass comparable to its stellar mass until $z = 9$. The virial mass accretion rate is generally consistent with the theoretical prediction for a $\Lambda$CDM cosmology in the Einstein–de Sitter (EdS) regime (Dekel et al. 2013):

where $s \simeq0.030\,\mathrm{Gyr}^{-1}$. The baryonic mass-flow history at $R_{\rm gal}$ is similarly smooth: the gas inflow rate gradually rises and remains larger than the outflow rate for most of the period shown, sustaining the gas reservoir within the galaxy. The SFR, averaged over both 5 Myr and 50 Myr, increases steadily with time. By $z = 9$, the SFR reaches $\sim50\,M_\odot\,\mathrm{yr}^{-1}$, comparable to both $\dot{M}_{\rm in}$ and $\dot{M}_{\rm out}$.

The projected UV half-light radius, $R_{\rm e}$ (magenta), grows from compact values at $z \gtrsim12$ to $\sim0.2$ kpc by $z = 9$, corresponding to 0.01–0.02 $R_{\rm v}$. The projected stellar half-mass radius, $R_{1/2,2\mathrm{D}}$ (grey), measured with the same source definition and estimator (§ 2.5), tracks $R_{\rm e}$ closely and is only marginally smaller at most epochs. The UV light and the stellar mass of the source are therefore distributed over similar scales, and the small residual offset is consistent with recent star formation extending slightly beyond the region containing half of the projected stellar mass. For comparison, Somerville et al. (2018) adopt projected mass-to-light size ratios of 0.83 and 0.87 for low-redshift optical disks and spheroids, respectively, similar to the ratio seen here. These optical values provide context rather than a conversion applied to our UV measurements, whose mass-to-light size ratio can change during individual star-formation episodes.

The 3D stellar half-mass radius, $R_{1/2}$ (black), is generally larger than its projected counterpart because projection reduces apparent stellar distances (Somerville et al. 2018; Behroozi et al. 2022). The shaded ranges in Figure 2 show the viewing-angle dependence. The offsets among the three radii reflect stellar-population and geometric effects, as well as the aperture and source definitions described in § 2.5. The 3D radius includes all stars within the spherical $R_{\rm gal}$ aperture, whereas image segmentation identifies the central source and can separate disconnected companions in the projected measurements. During the major galaxy merger at $z \simeq13.5$, the separation between the merging components produces a pronounced peak in $R_{1/2}$, while the projected radii remain more closely associated with the central component. This peak in $R_{1/2}$ traces the spatial extent of the combined stellar system during the interaction, rather than the size of an individual galaxy. Following the three radii together therefore helps distinguish the growth of an individual galaxy from the temporary contribution of companions within an aperture. The thick dotted grey line in the lower-right panel marks half the JWST/NIRCam F200W PSF FWHM. The projected UV radius becomes comparable to or larger than this reference scale near $z \sim9$, although its observational measurability also depends on signal-to-noise, PSF knowledge, and the assumed light profile.

Figure 3 presents the contrasting case of MAGE04. Unlike MAGE07, this system experiences a strongly merger-driven assembly history, as indicated by the lollipops in the bottom panels. Two major halo mergers occur around $z \sim12$, followed by two major galaxy mergers at $z \sim10.6$, which are visible in the top-middle panel as multiple bright components within the galactic region. Subsequently, the stellar mass increases rapidly, and the gas mass drops sharply, indicating efficient star formation triggered by wet mergers. By $z = 9$, the stellar mass reaches $\simeq0.31\,f_{\rm b}M_{\rm v}$, compared with $\simeq0.22\,f_{\rm b}M_{\rm v}$ for MAGE07.

Image not extracted

The virial mass accretion rate is broadly consistent with the theoretical expectation for a median halo mass accretion history (shown by the grey line in the middle panel) prior to the final major halo merger at $z \simeq11.7$. The corresponding gas inflow, outflow, and star-formation rates rise steadily, reaching $\sim100\,M_\odot\,\mathrm{yr}^{-1}$ during the major galaxy mergers at $z \simeq10.6$. Afterwards, the rates fluctuate around approximately constant levels, so that the final mass flow rates at $z = 9$ are comparable to those of MAGE07. The size evolution follows a similar pattern: $R_{\rm e}$ and $R_{1/2,2\mathrm{D}}$, which again remain close to each other, increase steadily before the major merger phase, with $R_{\rm e}$ exceeding 0.1 kpc already by $z \simeq11$. Subsequently, the galaxy size plateaus or undergoes net contraction, reflecting a wet compaction episode triggered by mergers that results in a final size at $z \sim9$ comparable to that of MAGE07. The transient peaks in $R_{1/2}$ again reflect multiple stellar components within $R_{\rm gal}$. Criterion (v) in § 2.2 reduces the contribution of such configurations to the statistical sample.

These two systems show that comparable sizes at $z \simeq9$ can arise from different assembly histories. MAGE07 grows relatively smoothly, whereas MAGE04 undergoes rapid merger-driven growth followed by a phase of gradual size contraction. Both end with $R_{\rm e} \sim0.2$ kpc, $R_{1/2} \sim0.2$–0.3 kpc, and comparable inflow, outflow, and star-formation rates. In both systems, the inflow and outflow rates at $R_{\rm gal}$ rise more steeply than the halo accretion rate, so the galaxy processes a growing fraction of the cosmological baryon supply, $f_{\rm b}M_{\rm v}$, without exceeding it. This suggests either more efficient penetration of fresh gas or the return through $R_{\rm gal}$ of gas ejected earlier. Meanwhile, the outflow rate falls from exceeding the SFR at early times to comparable values by $z = 9$, so the declining outflow-to-SFR ratio indicates that star formation becomes increasingly important relative to outward gas transport. This convergence motivates examining how the internal potential, feedback, and spatial distribution of star formation regulate sizes across the full sample.

The size–mass relation across redshift

We now examine the size–mass relation of galaxies at cosmic dawn using our statistical sample. Figure 4 shows the median 3D half-mass radius, $R_{1/2}$ (left panel), the 2D projected UV half-light radius, $R_{\rm e}$ (middle panel), and the corresponding ratio $R_{1/2}/R_{\rm v}$ (right panel) as a function of stellar mass in three redshift bins, with error bars spanning the 16th–84th percentiles. At a fixed stellar mass, we find little variation in $R_{1/2}$ and $R_{\rm e}$ across the three redshift bins, although the higher-redshift bins contain fewer data points at high stellar masses. The ratio $R_{1/2}/R_{\rm v}$ is slightly elevated at higher redshifts, primarily driven by a correspondingly smaller $R_{\rm v}$. These redshift dependencies are weaker than the mass dependence.

Image not extracted

We find a clear positive correlation between $R_{1/2}$ (as well as $R_{\rm e}$ and $R_{1/2}/R_{\rm v}$) and $M_\star$ for $M_\star< 5 \times10^{9}\,M_\odot$. At $M_\star\simeq5 \times10^{9}\,M_\odot$, the sizes reach $R_{1/2} \simeq300$ pc, $R_{\rm e} \simeq200$ pc, and $R_{1/2}/R_{\rm v} \simeq0.02$. Above this mass, the half-mass and UV half-light size–mass rela­tion tions become approximately flat or even turn over. This turnover is predominantly driven by the central galaxies in the four most massive targets (MAGE01–MAGE04), as other galaxies in the sample do not reach $M_\star> 5 \times10^9\,M_\odot$ by $z \simeq9$. The system MAGE04, shown in Figure 3, is one such example. While $R_{1/2}$ and $R_{\mathrm e}$ subsequently decrease by at most a factor of 2, the drop in $R_{1/2}/R_{\mathrm v}$ can exceed a factor of 2 due to the concurrent increase in $R_{\mathrm v}$ with stellar mass. At high mass, typical half-mass and UV half-light radii are 200–400 pc and 100–300 pc, respectively, for luminous objects at cosmic dawn.

To quantify the relation, we fit power laws of the form

directly to the galaxy–snapshot measurements (an unbinned fit) over the adopted stellar-mass range $\log_{10}(M_\star/M_\odot) \ge7.5$ for both $R_{1/2}$ and $R_{\mathrm e}$. The best-fit parameters are listed in Table 2. A single power law does not capture the high-mass behaviour. We therefore report separate fits over the intervals $7.5 \le\log_{10}(M_\star/M_\odot) < 9.25$ and $9.25 \le\log_{10}(M_\star/M_\odot) \le10.5$. These are independent fits over chosen mass intervals, rather than a fitted continuous broken power law.

Sample	$\log_{10}(M_\star/M_\odot)$	$N_{\mathrm{meas}}$	$\alpha$	$R_0$ [pc]	$\sigma_{\log R}$ [dex]
3D Half-Mass Radius ($R_{1/2}$)
$9 \le z < 11$, single	7.5–10.5	1585	$0.54 \pm0.01$	$127 \pm4$	0.28
$9 \le z < 11$, lower-mass interval	7.5–9.25	1520	$0.62 \pm0.02$	$154 \pm8$	0.28
$9 \le z < 11$, higher-mass interval	9.25–10.5	65	$-0.12 \pm0.04$	$340 \pm26$	0.10
$11 \le z < 13$	7.5–9.25	615	$0.66 \pm0.04$	$183 \pm18$	0.28
2D UV Half-Light Radius ($R_{\mathrm{e}}$)
$9 \le z < 11$, single	7.5–10.5	1585	$0.54 \pm0.01$	$79 \pm2$	0.22
$9 \le z < 11$, lower-mass interval	7.5–9.25	1520	$0.55 \pm0.02$	$81 \pm3$	0.22
$9 \le z < 11$, higher-mass interval	9.25–10.5	65	$0.04 \pm0.07$	$198 \pm25$	0.17
$11 \le z < 13$	7.5–9.25	615	$0.59 \pm0.03$	$97 \pm8$	0.22

In the left panel, we also compare our 3D half-mass radius–mass relation to those from THESAN-ZOOM (McClymont et al. 2025), THESAN-1 (Shen et al. 2024), FirstLight (Ceverino et al. 2026), and Lumina (Zier et al. 2026, Liu et al., to be submitted). Our galaxies are systematically smaller, and the offset grows towards low stellar mass. At $M_\star\simeq10^{7.75}\,M_\odot$, our median $R_{1/2} \simeq20$ pc lies more than an order of magnitude below the 0.5–1.5 kpc found in the four comparison suites, whereas by $M_\star\simeq10^{9-10.5}\,M_\odot$ our medians agree with FirstLight and lie within a factor of 2 of Lumina. At $M_\star\lesssim10^{7.5}\,M_\odot$, our half-mass radii are on the order of the spatial resolution of $\simeq10$ pc, which motivates the stellar-mass cut (criterion 1 in § 2.2) that selects only $M_\star\gtrsim10^{7.5}\,M_\odot$. The nature of these compact systems and the possible origins of the offset are discussed in § 6.5.1. Furthermore, our positive correlation at the low-mass end broadly aligns with THESAN-ZOOM, contrasting with the negative trend seen in THESAN-1, FirstLight, and Lumina over $10^{7.75} \lesssim M_\star/M_\odot\lesssim10^{9.25}$. In our sample, the half-mass relation turns over near $M_\star\simeq5 \times10^9\,M_\odot$, suggesting that the characteristic transition scale varies strongly between simulations (§ 6.5.2). We explore the implications of both this transition scale and our sub-grid feedback choices further in § 6.

In the middle panel, we compare our intrinsic 2D projected UV half-light radii with the corresponding dust-free values from THESAN-ZOOM, THESAN-1, FirstLight, and Lumina. Dust-attenuated sizes are not shown because they would not be comparable to our measurements. Even without dust, the comparison simulations lie within or above the range of JWST measurements (indicated by the navy and grey bands from Yang et al. 2025 and M24, respectively), whereas our sizes fall below it for $M_\star\lesssim10^{9.5}\,M_\odot$, with the same growing offset towards low mass as in the half-mass radii. The observed population is not confined to these bands, however.

Table 2. Power-law fits, Eq. (4), to the unbinned size–mass relation of central galaxies for both the 3D half-mass radius ($R_{1/2}$) and 2D UV half-light radius ($R_{\mathrm{e}}$), restricted to the adopted stellar-mass range $\log_{10}(M_\star/M_\odot) \ge7.5$. The quoted parameter uncertainties are formal standard errors from unweighted linear regression on galaxy–snapshot measurements; they do not include temporal correlations or uncertainty from the finite number of independent haloes. The residual scatter is given by $\sigma_{\log R}$, and $N_{\mathrm{meas}}$ counts galaxy–snapshot measurements.

M24 identify a compact subpopulation of 44 galaxies (≃ 13 per cent of their sample) that are marginally resolved or lie more than $1.5\sigma$ below their size–mass relation, with star-formation-rate surface densities of $\Sigma_{\mathrm{SFR}} > 10\,M_\odot\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$. The ten such sources at $z > 9$ (diamonds, seven with spectroscopic redshifts) have $R_{\mathrm{e}} \simeq10$–130 pc at $M_\star\simeq10^{7.4}$–$10^9\,M_\odot$ and fall between our medians and the mean relation.¹ About half of them, including both at $M_\star\gtrsim10^{8.5}\,M_\odot$, lie within or at the edge of our 16th–84th percentile range, and the most compact, with $R_{\mathrm{e}} \simeq10$–25 pc, coincide with our medians at $M_\star\simeq10^{7.6}$–$10^8\,M_\odot$, where the mean relation lies an order of magnitude above them. Our intrinsic sizes therefore overlap the compact tail of the observed distribution rather than its mean relation. The effects of dust attenuation and PSF on these comparisons, which can be substantial for compact galaxies, are discussed in § 6.6.1 and § 6.8.

The redshift evolution is further detailed in Figure 5, where $R_{1/2}$ (left panel), $R_{\mathrm{e}}$ (middle panel), and $R_{1/2}/R_{\mathrm{v}}$ (right panel) are shown as a function of $z$ across different stellar-mass bins. At a fixed stellar mass, $R_{1/2}$, $R_{\mathrm{e}}$, and $R_{1/2}/R_{\mathrm{v}}$ exhibit minimal evolution with cosmic time, reflecting the weak redshift dependence observed in Figure 4. This stability broadly aligns with other simulations (THESAN-ZOOM, FIRSTLIGHT, LUMINA) and JWST constraints (M24), though our intrinsic galaxy sizes remain systematically more compact. The offset is again largest in the lowest-mass bin: for $7.5 \le\log_{10}(M_\star/M_\odot) < 8.5$, FIRSTLIGHT and LUMINA maintain $R_{1/2} \sim0.5$–$0.7$ kpc over $9 \lesssim z \lesssim12$, compared with $\sim25$ pc in MAGE, whereas for $8.5 \le\log_{10}(M_\star/M_\odot) < 9.5$ our half-mass radii @page 10/30 starts_mid=0 ends_mid=1 approach those of FirstLight at $z \gtrsim10$. The compact sources of M24 at $z > 9$ (diamonds in the middle panel) have $R_{\rm e} \simeq10$–$130$ pc over $z \simeq9$–13, within the range spanned by our two lower-mass bins at the same redshifts and well below the band for their full sample at $M_\star= 10^8 M_\odot$. The medians in the two higher stellar-mass bins indicate a typical size ratio of $R_{1/2}/R_{\rm v} \simeq0.015$ for galaxies with $M_\star> 10^{8.5} M_\odot$, comparable to $R_{\rm e}/R_{\rm v} \simeq0.015$ of Ono et al. (2025) at $z \sim10$–16, shown as the grey dashed line in the right panels of Figure 4 and Figure 5, as well as the nearly constant half-light ratios of $R_{\rm e}/R_{\rm v} \simeq0.01$–0.035 measured from $z \sim0$ to $z \sim10$ (Shibuya et al. 2015; Yang et al. 2025) and $\simeq0.015$–0.02 inferred at $z \lesssim3$ (Kravtsov 2013; Somerville et al. 2018), despite using different size definitions.

Image not extracted

The principal result is a positive low-mass size–mass relation with weak redshift evolution, followed by a flattening or turnover in the half-mass and UV half-light size–mass relations at high mass. The transition occurs at a higher stellar mass than in several comparison simulations, while our intrinsic sizes have a lower normalization, most strongly at $M_\star\lesssim10^8 M_\odot$, where our galaxies are compact and star-cluster-like rather than extended. Relative to observations, our intrinsic UV sizes lie below the mean size–mass relation but overlap the compact subpopulation of M24.

The size–magnitude relation

Because stellar-mass estimates at $z \geq10$ are uncertain, we now examine the size–magnitude relation at $9 \leq z < 15$, which can be compared more directly with observations. Figure 6 shows the median 2D projected UV half-light radius, $R_{\rm e}$, as a function of rest-frame UV magnitude for all central galaxies across three redshift bins, compared with the spectroscopically confirmed $z \geq10$ sample of Roberts-Borsani et al. (2026), split by whether C IV emission is detected in their NIRSpec prism spectra, and with the empirical relation of M24. Absolute UV magnitudes in MAGE are evaluated at $\lambda= 1500$ Å from the selected source, using the aperture conventions described in § 2.5. These monochromatic magnitudes differ slightly from the F200W band used for the sizes, which samples a redshift-dependent rest-frame wavelength, but we expect this to have

Image not extracted

a negligible effect. They are intrinsic, with no dust attenuation, PSF effects, or noise included.

At fixed absolute UV magnitude, galaxies at higher redshift are more compact, a stronger redshift evolution than seen in the size–mass relation. This is consistent with the younger stellar populations at higher redshift are brighter, so a given $M_{\mathrm{UV}}$ corresponds to a lower stellar mass. In each redshift interval, brighter galaxies are generally larger. The slope is steep between $M_{\mathrm{UV} }=-19$ and $-21$, then flattens at the bright end, where it is more consistent with the slope of the typical sample of M24, although our galaxies lie close to the lower envelope of the distribution. At the faint end ($M_{\mathrm{UV}}>-20$), our galaxies fall below the scatter of the typical sample and populate the compact tail of the observed size distribution, coinciding with the population of C iv-strong sources (Roberts-Borsani et al. 2026).

Part of the offset is definitional: M24 report major-axis effective radii, whereas both Roberts-Borsani et al. (2026) and we use circularized half-light radii. For an axis ratio $q$, a major-axis radius exceeds the corresponding circularized radius by a factor of approximately $q^{-1/2}$, or $\gtrsim0.1$ dex for $q\lesssim0.6$. This convention alone therefore shifts their relation upward relative to ours, although a precise correction would require matched radius estimators. Dust attenuation, the finite angular resolution of the observations, our limited sample size, and our exclusion of interacting systems also affect this comparison. Their effects on the observed sizes and magnitudes are discussed in § 6.6.1. Finally, for galaxies fainter than $M_{\mathrm{UV}}\simeq-19$, the measured sizes approach the numerical resolution limit of the simulations.

In-situ vs. ex-situ growth: how did early galaxies assemble their stars?

In § 4 and § 5, we explore the origin of size evolution at cosmic dawn. We follow the central galaxies of the 11 target haloes, whose main-progenitor histories allow a more detailed analysis than the full statistical sample. Motivated by the turnover in galaxy size at $M_{\star}\simeq5\times10^{9}\,M_{\odot}$ found in § 3.2, and by a second characteristic mass at $M_{\star}\simeq5\times10^{8}\,M_{\odot}$ that emerges from the kinematic analysis in § 5, we divide our sample of 11 central galaxies into three groups based on their final stellar mass at $z=9$. In order of increasing mass, these are a low-mass group ($M_{\star}\leq5\times10^{8}\,M_{\odot}$; MAGE09–11), an intermediate-mass group ($5\times10^{8}\,M_{\odot}<M_{\star}\leq5\times10^{9}\,M_{\odot}$; MAGE05–08), and a high-mass group ($M_{\star}>5\times10^{9}\,M_{\odot}$; MAGE01–04). This grouping follows each system according to its final mass. An individual progenitor can cross several instantaneous mass intervals during its evolution. We relate these intervals to the structural and kinematic transitions in § 5.

We compare the spatial distributions of in-situ and ex-situ stars defined in § 2.2 to assess their contributions to size evolution. Figure 7 provides a visual overview at $z=9$, showing face-on stellar surface-density maps of the 11 target central galaxies within $\pm4R_{1/2}$ and contours of the projected in-situ fraction at 0.5, 0.7, and 0.9. In ten of the eleven galaxies, the 0.7 contour encloses the central stellar body, and in several systems (e.g. MAGE01, MAGE07, MAGE10, and MAGE11) the in-situ fraction exceeds 0.9 over much of the inner region. In-situ stars therefore dominate the central stellar mass in ten of the eleven systems, spanning a broad range of masses and merger histories. The predominantly in-situ central mass of MAGE04, despite its recent major mergers, is consistent with the high gas fractions during the interaction (bottom-left panel of Figure 3): stars formed locally can contribute more central mass than the directly accreted stellar population.

Image not extracted

The one exception is MAGE08, which experienced a major merger within the last 50 Myr. Its in-situ fraction exceeds 0.5 only in small central patches, so directly accreted stars dominate its stellar body at $z=9$. In the other systems dominated by in-situ stars, ex-situ stars concentrate in the outskirts, where the 0.5 contour isolates satellites and tidal debris, and in localized central dips in the in-situ fraction such as the small closed contour at the centre of MAGE01. These dips likely trace the remnants of accreted companions that have sunk to the core via dynamical friction (Chandra­se­khar 1943; Barnes & Hernquist 1992; Boylan-Kolchin et al. 2008). The projected maps thus show predominantly in-situ growth in the central regions, with ex-situ stars dominating the stellar body only in the one galaxy observed immediately after a major merger.

The left panel of Figure 8 presents the radial profiles of the 3D in-situ stellar mass fraction at $z=9$ for all central galaxies, split into the three mass groups. In all three groups, the median in-situ fraction is $f_{\mathrm{in-situ}}\gtrsim0.7$ within $R_{1/2}$ and declines mildly towards the centre. This central dip is likely driven either by ex-situ stars that have sunk inward through dynamical friction, as suggested by the closed central contour of MAGE01 in Figure 7, or by more recent in-situ star formation near $R_{1/2}$, as shown in the right panel of Figure 8. Beyond $\simeq1.5$–$2R_{1/2}$, the median of the low-mass group remains at $\simeq0.8$ out to $5R_{1/2}$, albeit with only three galaxies and a large spread among them, whereas the in-situ fraction of the intermediate- and high-mass groups drops rapidly, reaching $f_{\mathrm{in-situ}}\simeq0.2$–$0.4$ at $R\simeq5R_{1/2}$.

Image not extracted

The galaxy-to-galaxy scatter within $R_{1/2}$ also differs between the groups: the three low-mass galaxies span $f_{\mathrm{in-situ}}\simeq0.7$–1. The intermediate-mass group includes one galaxy (MAGE08) with $f_{\mathrm{in-situ}}\simeq0.3$–0.5 throughout its stellar body due to a very recent major merger within 50 Myr, and the high-mass group exhibits the smallest scatter, with all four galaxies confined to $f_{\mathrm{in-situ}}\simeq0.8$–1. The low scatter within the high-mass group is consistent with a rapid central stellar mass buildup (discussed below), which anchors their inner $f_{\mathrm{in-situ}}$ near unity despite their different assembly histories in this sample. Finally, the increasing scatter at large radii reflects stochastic contributions from mergers, accreted satellites, and tidal debris, which predominantly populate the galactic outskirts.

The fact that in-situ stars dominate the central stellar mass does not by itself establish that they control size evolution. To investigate this, we compare the half-mass radii of the total, in-situ, and ex-situ stellar populations in Figure 9, showing the three mass groups in separate panels. In all three groups, the total $R_{1/2}$ broadly follows the in-situ radius and is smaller than the ex-situ radius. The in-situ component exhibits a smoother growth history and closely tracks the total $R_{1/2}$ for the majority of the evolution. In the high-mass group, the two become nearly indistinguishable after $t\simeq440$ Myr. In contrast, the ex-situ $R_{1/2}$ is a factor of $\sim2$–4 larger and highly stochastic, driving several transient peaks in the overall size (e.g., at $t\simeq525$ Myr in the low-mass group and at $t\simeq320$ Myr in the intermediate-mass group). This behaviour supports a scenario where in-situ star formation governs the baseline structural evolution. Ex-situ stars, accreted via mergers, settle into a more extended distribution and their smaller mass contribution generally limits their effect on the total $R_{1/2}$.

Image not extracted

The growth histories differ markedly between the groups. The low-mass group grows only slowly, from $\simeq15$ pc to $\simeq50$ pc, sizes that are only a few times the $\sim10$ pc spatial resolution. The intermediate-mass group exhibits the steepest temporal evolution, growing steadily by an order of magnitude, from $\simeq30$ pc at $t\simeq300$ Myr to $\simeq250$ pc at $z=9$. In the high-mass group, the in-situ $R_{1/2}$ grows from $\simeq80$ pc at $t\simeq300$ Myr to $\simeq250$ pc by $t\simeq440$ Myr and flattens thereafter, so that the intermediate-mass group catches up and reaches comparable sizes by $z = 9$. These trends align with the size–mass relation presented in Figure 4, which exhibits a steep slope at the low-mass end that flattens—and even turns over—above $M_\star\simeq5 \times10^9 M_\odot$.

The radial ages of in-situ stars provide a complementary view of their formation history. The right panel of Figure 8 shows their radial age profiles at $z = 9$. All three groups exhibit an inside-out formation signature, with the oldest stars at the centre and the youngest at $R \gtrsim R_{1/2}$, but the shape of the age profile changes systematically with mass. In the low-mass group, the median age is nearly flat at $\simeq90$–110 Myr throughout $R < R_{1/2}$ and then drops abruptly to $\simeq50$–65 Myr outside $R_{1/2}$. The intermediate-mass group shows a smooth, shallow age gradient and has the lowest median stellar ages of the three groups at all radii within $R_{1/2}$. Its median age declines from $\simeq75$ Myr at the centre to $\simeq45$ Myr at 1–2 $R_{1/2}$, consistent with the steady, ongoing size growth of this group in Figure 9. The high-mass group shows a smooth but steeper gradient. The median stellar age in its innermost region is 110 Myr, corresponding to an epoch of $t \simeq440$ Myr ($z \simeq10.6$). This roughly coincides with the epoch when MAGE04 experiences a merger-driven peak in star formation (Figure 3). Moving outward, the median stellar age declines smoothly, reaches a minimum of $\simeq50$ Myr near $R_{1/2}$, and remains low out to $5R_{1/2}$, indicating a more recently formed population extended across the radial region. The older centre and younger outskirts of the high-mass group are consistent with more pronounced inside-out quenching following the rapid central stellar buildup: star formation declines in the compact core while continuing at larger radii. These galaxies, with $M_\star> 5 \times10^9 M_\odot$, have also undergone the size turnover and developed substantial rotational support (§ 5).

The birth radii of the in-situ stars present at $z = 9$ provide a retrospective measure of their formation sites (Figure 10). This analysis is restricted to stars retained within the selected galaxy at that epoch. We bin these stars by the stellar mass of the main progenitor at their formation time, $M_{\star}(t_{\mathrm{birth}})$ (§ 2.2). Since the 11 galaxies reach a given mass at different epochs, this places them on a common evolutionary sequence and allows a direct comparison with the size–mass relation (Figure 4) and with the transition masses identified in § 5. Panel (a) displays the absolute birth radii. The median $R_{\mathrm{birth}}$ is $\simeq100$–$200$ pc for $M_{\star} \lesssim10^{9}\,\mathrm{M}_{\odot}$, rises to $\simeq400$–$500$ pc at $M_{\star} \simeq(3$–$5) \times10^{9}\,\mathrm{M}_{\odot}$, and then flattens, and even declines slightly, at higher masses, mirroring the turnover of the size–mass relation. Panel (b) compares the birth radii to the half-mass radii at the time of formation. At $M_{\star} \lesssim5 \times10^{8}\,\mathrm{M}_{\odot}$, the median birth radius exceeds $R_{1/2}$ by a factor of $\sim2$, with a galaxy-to-galaxy range spanning more than an order of magnitude: stars form predominantly beyond the existing stellar half-mass radius, and this geometric offset favours expansion of $R_{1/2}$ and the steep slope of the size–mass relation at the low-mass end observed in Figure 4. Above $M_{\star} \simeq5 \times10^{8}\,\mathrm{M}_{\odot}$, the ratio converges to $R_{\mathrm{birth}}/R_{1/2} \simeq1$–$1.5$. The galaxy-to-galaxy range narrows markedly, indicating that star formation tracks the half-mass radius more closely, although this convergence precedes the turnover near $5 \times10^{9}\,\mathrm{M}_{\odot}$.

Image not extracted

We explore stellar redistribution by comparing these birth radii to the galactocentric radii of the same selected stars at $z = 9$ (panel (c) of Figure 10). The median $R_{\mathrm{birth}}/R_{\mathrm{now}}$ is generally $\simeq1$–$1.5$. It peaks at $\sim2$ for stars formed at $M_{\star} \simeq5 \times10^{9}\,\mathrm{M}_{\odot}$ and returns to unity for those formed at $M_{\star} \gtrsim10^{10}\,\mathrm{M}_{\odot}$. The galaxy-to-galaxy range is large for stars formed at $M_{\star} \lesssim10^{8}\,\mathrm{M}_{\odot}$, consistent with the more chaotic assembly of the low-mass progenitors, but narrows considerably at higher masses. These ratios indicate modest net decreases in galactocentric radius for most cohorts, with larger decreases near the turnover. Finally, the normalized birth radius, $R_{\mathrm{birth}}/R_{\mathrm{v}}$ (panel (d)), fluctuates around $\simeq0.02$ without a systematic trend up to $M_\star\simeq5 \times10^9\,M_\odot$, and then declines to $\simeq0.01$ at $M_\star\gtrsim10^{10}\,M_\odot$.

Overall, a coherent picture emerges along the stellar-mass sequence. Below $M_\star\simeq5 \times10^8\,M_\odot$, the star-forming radius inferred from stellar birth sites tracks the halo at $\simeq0.02\,R_{\mathrm{v}}$ and lies well outside $R_{1/2}$, driving rapid expansion of the half-mass radius toward this outer boundary. Between $M_\star\simeq5 \times10^8\,M_\odot$ and $\simeq5 \times10^9\,M_\odot$, the two radii have converged and grow together. Above $M_\star\simeq5 \times10^9\,M_\odot$, the star-forming radius saturates around $\sim400$ pc while the halo continues to grow, so that $R_{\mathrm{birth}}/R_{\mathrm{v}}$, and hence $R_{1/2}/R_{\mathrm{v}}$, decline. This coincides with the turnover of the size–mass relation in Figure 4. This sequence suggests a transition from size growth that tracks the halo to stronger internal regulation at high stellar mass. We return to the physical origin of these transition masses in § 6.

Disk Formation in Massive Galaxies

We now connect the kinematics of the star-forming gas, the newly formed stars, and the full stellar population to size evolution, using the rotation-to-dispersion ratio $v/\sigma$ defined in § 2.3. In the diagnostics below, we measure the inflow rate across $R_{\mathrm{gal}}$, the outflow rates across $R_{1/2}$, $R_{\mathrm{gal}}$, and $R_{\mathrm{v}}$, and the SFRs within $R_{1/2}$ and $R_{\mathrm{gal}}$. Figure 11 presents the time evolution of $v/\sigma$ for the star-forming gas, for stars younger than 10 Myr, and for the full stellar population, divided into the three final-mass groups defined in § 4. The low-mass group fluctuates strongly without a sustained increase. Apart from transiently high values at $t \lesssim330\,\mathrm{Myr}$ ($z \gtrsim13$), the median $v/\sigma$ of its young stars varies between $\simeq0.6$ and 1.8, with a temporary peak near $t \simeq500\,\mathrm{Myr}$ that has fallen back to $\simeq1.1$ by $z = 9$, while the star-forming gas and the full stellar population remain at $v/\sigma\lesssim1$. High rotational support for young stars can therefore occur transiently in low-mass systems and does not by itself establish a lasting disk. The intermediate-mass group evolves more gradually: the median $v/\sigma$ of its young stars remains near $\simeq0.8$–1.4 for most of the interval and reaches $\simeq1.5$ by $z = 9$, the star-forming gas stays at $\simeq0.5$–0.9 apart from isolated excursions and reaches $\simeq1$–1.2 only after $t \simeq520\,\mathrm{Myr}$, and the full stellar population remains at $\simeq0.5$–0.8. Only the high-mass group exhibits a sharp, sustained kinematic transition near $t \simeq470\,\mathrm{Myr}$ ($z \simeq10$). Before the transition, all three tracers remain below unity, mostly at $v/\sigma\simeq0.2$–0.8, for $\simeq150\,\mathrm{Myr}$. Within a few tens of Myr thereafter, the median $v/\sigma$ of the young stars rises from $\lesssim1$ to $\gtrsim2$ and that of the star-forming gas from $\simeq1$ to $\simeq1.6$. By $z = 9$, the young stars reach $v/\sigma\simeq3$ and the star-forming gas $\simeq1.6$, whereas the full stellar population responds more slowly, rising from $\simeq0.4$ at the transition to $\simeq1$ at $z = 9$. Rotational support is thus established first in the gas and the newly formed stars and spreads only gradually to the stellar body as a whole. The transition occurs approximately 30 Myr after the formation epoch implied by the median central stellar age, $t \simeq440\,\mathrm{Myr}$ (Figure 8, right panel), consistent with central stellar buildup preceding rapid rotational support.

Image not extracted

Panel (d) of Figure 12 shows the mass dependence of $v/\sigma$ for the star-forming gas, for stars younger than 10 Myr, and for the full stellar population. The star-forming gas and the young stars track each other below the transition, with median $v/\sigma\simeq0.5$–$1$, and both rise above unity near $M_\star\simeq5 \times10^{9}\,M_\odot$. The young stars reach $v/\sigma\simeq2$–$3$, above the $\simeq1.5$ of the gas. Part of this difference reflects the thermal and turbulent support included in $\sigma_{\mathrm{SFgas}}$, which the stars do not inherit. The median full-population ratio remains approximately at or below unity, with $v/\sigma\lesssim1$ even at the highest masses, so the rotational support is confined to the newly formed component inside a stellar body that stays dispersion-dominated. Above $M_\star\simeq1.5 \times10^{10}\,M_\odot$, only MAGE01 remains in the sample. It experiences a sharp decline in rotational support due to mergers, followed by a recovery, eventually reaching its peak $(v/\sigma)_{\star,10} \gtrsim3$ at $M_\star> 3 \times10^{10}\,M_\odot$. The rapid onset of rotational support near $M_\star\simeq5 \times10^{9}\,M_\odot$ is therefore common to the gas and the young stars.

Image not extracted

The stellar-mass threshold for disk formation coincides with the turnover in the half-mass radius (panel (c)). Although this panel isolates the 11 central galaxies without binning by redshift, the weak redshift dependence permits comparison with the broader morphological trends shown in the left panel of Figure 4. Panel (b) shows the SFR-weighted rotation speed and velocity dispersion of the star-forming gas, together with the circular velocity at the half-mass radius and the halo virial velocity. The effective dispersion tracks the halo velocity scale over the full mass range, $\sigma_{\mathrm{SFgas}} \simeq0.6\,v_{\mathrm{v}}$, with no systematic change across the transition. The rotation speed behaves differently. Below the transition, $v_{\mathrm{rot,SFgas}} \simeq0.3$–$0.5\,v_{\mathrm{c,tot}}(R_{1/2})$ and lies below $\sigma_{\mathrm{SFgas}}$. Near $M_\star\simeq5 \times10^{9}\,M_\odot$ where $v_{\mathrm{c,tot}}$ exceeds $v_{\mathrm{v}}$, it rises steeply to $\simeq0.6$–$0.7\,v_{\mathrm{c,tot}}(R_{1/2})$ and overtakes $\sigma_{\mathrm{SFgas}}$. The rise in $v/\sigma$ at the transition is therefore driven by the emergence of coherent rotation at a dispersion that continues to follow $v_{\mathrm{v}}$, rather than by a drop in dispersion. Even after the transition the gas does not rotate at the circular speed at $R_{1/2}$. $v_{\mathrm{rot,SFgas}}$ is an SFR-weighted mean over the gas within $R_{\mathrm{gal}}$, which includes non-circular and misaligned motions and gas at radii where the circular speed differs from its value at $R_{1/2}$, so the ratio should be read as a measure of coherent rotation rather than of the departure from circular orbits at a fixed radius.

To understand why $v_{\mathrm{c,tot}}$ eventually exceeds $v_{\mathrm{v}}$, we examine the stellar, dark matter, and gas contributions in panel (a). This steep rise in $v_{\mathrm{c,tot}}$ is driven by a concurrent surge in the stellar component, $v_{\mathrm{c,star}}$ (orange line), which exceeds $v_{\mathrm{c,dm}}$ from $M_\star\simeq5 \times10^{8}\,M_\odot$ onward, signalling a phase of rapid central mass buildup prior to disk formation. We adopt a characteristic velocity $v_{\mathrm{SN}}$, the critical virial velocity of Dekel & Silk (1986) shown by the horizontal dotted line. It is defined by equating the energy deposited in the interstellar gas by SNe during the first burst of star formation, after accounting for radiative losses from the remnants, with the gravitational binding energy of the gas in the halo. Haloes with $v_{\mathrm{v}} < v_{\mathrm{SN}}$ can lose a substantial fraction of their gas to SN-driven winds, whereas deeper potential wells retain it. Because the deposited SN energy scales with the stellar mass formed, nearly independently of the gas density, the resulting threshold depends only weakly on the collapse epoch. We use $v_{\mathrm{SN}}$ as a characteristic potential-depth scale associated with effective SN feedback, not as a measured wind launch speed or an exact escape criterion.

We find that the virial velocity $v_{\mathrm{v}}$ and the central circular velocity $v_{\mathrm{c,tot}}$ surpass $v_{\mathrm{SN}}$ at $M_\star\simeq10^{8}\,M_\odot$ and $M_\star\simeq5 \times10^{8}\,M_\odot$, respectively. Near $M_\star\simeq10^{8}\,M_\odot$, the half-mass radius begins to grow more rapidly as the halo virial speed exceeds $\simeq100\,\mathrm{km\,s^{-1}}$. This is consistent with a scenario in which gas is ejected to the outer circumgalactic medium (CGM), mixes with high-angular-momentum gas, and falls back onto the outskirts of the galaxy. In this case, the angular momentum of the CGM gas may help the forwards track halo growth, so the star-forming radius scales with the halo radius (§ 4). Near $M_\star\simeq5 \times10^{8}\,M_\odot$, the total circular speed at $R_{1/2}$ exceeds $v_{\mathrm{SN}}$, consistent with stronger retention of SN-ejected gas by the central potential. In this interpretation, faster replenishment sustains central star formation and cumulative stellar buildup, while the gas fraction declines and the stellar contribution overtakes that of dark matter within $R_{1/2}$. Size growth slows as the star-forming radius converges on the half-mass radius (Figure 10), until the final size turnover at $M_\star\simeq5 \times10^{9}\,M_\odot$.

Panel (e) supports this picture through the outflow rates measured at $R_{1/2}$, $R_{\rm gal}$, and $R_{\rm v}$, together with the inflow rate at $R_{\rm gal}$. We normalize these rates by the SFR within $R_{\rm gal}$, averaged over the past 50 Myr. The normalization expresses each rate relative to the galaxy-wide SFR and removes the cosmic increase in the rates as the galaxy grows. The outflow rates include all gas moving outward through each shell (§ 2.2), without tracing its origin or applying a velocity cut, so they include bound turbulent and orbital motions and, at the outer shells, gas associated with satellites and tidal debris that never belonged to the central galaxy. They are therefore crude measures of the competition between feedback and gravity, and rely on their ordering and trends rather than on their absolute values. The overall trend, highlighted by the arrow in panel (e), is towards less large-scale outflow at higher mass: the rate across $R_{\rm v}$ falls below that across $R_{\rm gal}$, which in turn falls below that across $R_{1/2}$, an ordering consistent with increasing confinement of outward-moving gas to smaller radii as the potential deepens. The mass loading across $R_{\rm v}$ declines from $\simeq1.5$–$2$ at $M_\star\lesssim5 \times10^{8}\,M_\odot$ to $\simeq0.4$ at $M_\star\gtrsim3 \times10^{9}\,M_\odot$, as expected if a deepening potential retains a growing share of the ejected gas, and lies a factor of $\sim2$ below that across $R_{\rm gal}$ over the full mass range.

The central star-formation fraction, defined as the SFR within $R_{1/2}$ divided by the SFR within $R_{\rm gal}$, both averaged over the past 5 Myr (dashed orange curve in panel (e)), stays at $\simeq0.4$, within $0.3$–$0.55$, up to and across the three transitions, and subsequently declines in the final bins at $M_\star\gtrsim2.5 \times10^{10}\,M_\odot$, where $v_{\rm c,gas}$ shows a similar sharp decline. The decline of the central star-formation fraction, the low central gas fraction, and the younger outskirts in the age profiles of Figure 8 together suggest the onset of inside-out quenching, although this decline is traced only by MAGE01. We assemble these trends into a sequence of regimes in § 6.7 (Figure 15).

Image not extracted

Discussion

Connection to wet compaction events

The central mass buildup, size turnover, disk formation, and inside-out quenching identified in § 4 and § 5 share several features with wet compaction at cosmic noon. In that picture, dissipative gas inflow builds a dense, star-forming core, often termed a “blue nugget”, and changes both the central mass composition and the galaxy kinematics (Dekel & Burkert 2014; Zolotov et al. 2015; Tacchella et al. 2016; Lapiner et al. 2023). The end of compaction and the onset of disk formation occur at a characteristic stellar mass of $\sim10^{10}\,\mathrm{M}_{\odot}$ (Lapiner et al. 2023). Blue dashed line in the bottom panel of Figure 1, only a factor of two above the characteristic stellar mass at which we find rapid disk formation.

In these studies, compaction is typically triggered by wet mergers, counter-rotating streams, or violent disk instability, though Lapiner et al. (2023) also identified gas recycling as a viable driver. Among the four massive galaxies in our sample that undergo compaction, MAGE04 is likely driven by wet mergers (Figure 3). In contrast, the remaining three systems (MAGE01–MAGE03) exhibit no evidence of prior major mergers, counter-rotating streams, or violent disk instability preceding disk formation. The compaction-like behavior in these MAGE galaxies is therefore likely driven by rapid internal gas recycling, aligning with the broader theoretical framework of compaction. Our interpretation emphasizes the internal conditions that confine outflows, permit rapid recycling, and allow central gas accumulation to persist long enough to become structurally significant, without attributing the transition to any particular inflow trigger. Stronger inflows induced by wet mergers, for instance, would accelerate central gas accumulation. The proposed sequence begins when $v_{\mathrm{c,tot}}(R_{1/2}) > v_{\mathrm{SN}}$, marking a shift from halo-scale to galaxy-scale recycling, and is followed by rapid growth of the central potential. Cosmological accretion continues to supply gas, and interactions may shape individual bursts, but the proposed regulator of the characteristic transition is the evolving balance between central gravity and feedback.

In more detail, once $v_{\mathrm{c,tot}}(R_{1/2})$ exceeds the characteristic feedback scale $v_{\mathrm{SN}}$ adopted in § 5, at $M_{\star} \simeq5 \times10^{8}\,\mathrm{M}_{\odot}$, feedback becomes less effective at removing gas from the galaxy: the outflow rate across $R_{\mathrm{gal}}$ becomes lower than the rate across $R_{1/2}$ (Figure 12). This ordering is consistent with outflows failing to escape the galaxy and recycling there, rather than reaching the CGM before returning to the outskirts, although shell fluxes alone do not establish those trajectories. The ratio of the outflow rate across $R_{1/2}$ to the estimated net inflow across $R_{\mathrm{gal}}$ rises from $\lesssim0.1$ at $M_{\star} \simeq10^{8}\,\mathrm{M}_{\odot}$ to approximately unity near $5 \times10^{9}\,\mathrm{M}_{\odot}$, suggesting that recycling may contribute as much as fresh accretion to the central gas supply. Note that bursty star formation in individual galaxies is not confined to this regime. Galaxies between $10^{8}\,\mathrm{M}_{\odot}$ and $5 \times10^{8}\,\mathrm{M}_{\odot}$ also undergo bursts, but the proposed scenario predicts that gas expelled from these galaxies returns on the longer halo-scale timescale, leading to longer intervals before the gas supply is restored and therefore a lower burst duty cycle. Recycling within the galaxy shortens the replenishment timescale, so the central gas supply is restored more quickly after each star-forming episode. The expected signature is a shorter burst cycle with a higher duty cycle: star formation recurs more frequently, and a gas supply maintained between bursts reduces the contrast between burst peaks and the inter-burst level without eliminating the central enhancements. A follow-up analysis of star-formation variability using MAGE indicates reduced burstiness above $M_{\star} \simeq5 \times10^{8}\,\mathrm{M}_{\odot}$ (Yao et al., in preparation), consistent with more effective gas confinement and shorter replenishment timescales. The accumulating stellar mass deepens the potential and further promotes gas retention, potentially sustaining efficient central star formation until star formation consumes the gas and/or redistribution reduces the central gas supply. A related consequence of this central buildup may be the suppression of global “breathing modes”, in which gas coherently expands and contracts across the galaxy. A more centrally concentrated mass profile causes radial oscillations at different radii to lose coherence, promoting gas mixing and facilitating disk formation (Hopkins et al. 2023).

Although individual galaxies undergo bursts in both regimes, the rapid stellar buildup near $M_{\star} \simeq5 \times10^{9}\,\mathrm{M}_{\odot}$ is the phase in our massive systems that coincides with both the size turnover and the sharp increase in rotational support. By this stage, the accumulated central mass has raised $v_{\mathrm{c,tot}}(R_{1/2})$ above $v_{\mathrm{v}}$. The condition for retaining gas within the galaxy is therefore met before the condition associated with rapid disk formation: the central stellar buildup between $5 \times10^{8}\,\mathrm{M}_{\odot}$ and $5 \times10^{9}\,\mathrm{M}_{\odot}$ assembles the core, and the rapid buildup near $5 \times10^{9}\,\mathrm{M}_{\odot}$ marks the point at which the accumulated mass changes both the size evolution and the kinematics. In this picture, the compaction-like phase in MAGE may begin near $M_{\star} \simeq5 \times10^{8}\,\mathrm{M}_{\odot}$, where the growth of $R_{1/2}$ slows from rapid expansion. The size turnover near $5 \times10^{9}\,\mathrm{M}_{\odot}$ may therefore reflect the continued decrease in the size-slope during central buildup, rather than a separate event. Note that the characteristic mass of $\sim10^{10}\,\mathrm{M}_{\odot}$ in Lapiner et al. (2023) marks the end of compaction, after which the galaxy’s half-mass radius begins to increase again. Our simulations have not yet reached this phase of renewed size growth (§ 6.7.2). The size turnover and the onset of rotational support in MAGE may thus be two consequences of the same central stellar buildup, although they need not occur at exactly the same stellar mass.

The role of gas recycling in this buildup remains an interpretation: our shell-flux measurements do not distinguish fresh accretion from returning gas or track its angular-momentum loss. Establishing the relative contributions of these processes requires a dedicated analysis of gas trajectories. The connection between central mass concentration and rotational support can nevertheless be examined directly, as we do next in § 6.2.

The conditions for disk formation

In § 5, the rapid emergence of rotational support in our massive galaxies coincides with $v_{\mathrm{c,tot}}(R_{1/2}) > v_{\mathrm{v}}$ at a stellar mass of approximately $5 \times10^{9}\,\mathrm{M}_{\odot}$. This comparison of the central circular speed with the halo velocity scale provides an indirect diagnostic of disk formation. The connection depends on the actual rotation speed and velocity dispersion of the gas and stars. Writing $v_{\mathrm{c}} = v_{\mathrm{c,tot}}(R_{1/2})$, the rotation-to-dispersion ratio can be written as the product

where $v_{\mathrm{rot}}$ denotes the mean rotation speed of the population whose dispersion is $\sigma$. For the star-forming gas, the relevant quantities are shown in panel (b) of Figure 12. The velocity dispersion remains at $\sigma_{\mathrm{SFGas}} \simeq0.6\,v_{\mathrm{v}}$ across the mass range, while $v_{\mathrm{rot,SFGas}}/v_{\mathrm{c}}$ rises from $\simeq0.4$ below the transition to $\simeq0.7$ above it. Once the disk has formed, $v/\sigma$ is therefore approximately $v_{\rm c}/v_{\rm v}$, to within $\sim20$ per cent. The substantially steeper rise in $v_{\rm rot,SFgas}$ relative to $v_{\rm c,tot}$ warrants further investigation. One plausible mechanism is the suppression of the breathing mode, which facilitates gas mixing so that $v_{\rm rot,SFgas}$ approaches $v_{\rm c,tot}$ (Hopkins et al. 2023). A definitive numerical verification is deferred to future work. Under the approximate assumptions $v_{\rm rot} \sim v_{\rm c}$ and $\sigma\sim v_{\rm v}$, the criterion for rotational support, $v/\sigma> 1$, becomes $v_{\rm c}/v_{\rm v} > 1$. This condition is met near the transition mass of $M_{\star} \simeq5 \times10^{9}\,M_{\odot}$. Since $v_{\rm c,tot}(R_{\rm v}) = v_{\rm v}$, this inequality implies a centrally peaked rotation curve that declines with radius. Because stars dominate the enclosed mass at these stellar masses, the more conservative diagnostic $v_{\rm c,star}(R_{1/2}) > v_{\rm v}$ can be used. Under spherical symmetry, this criterion becomes

This expression relates the stellar-to-halo mass ratio to the stellar-to-halo radius ratio, and identifies compact galaxies with a large stellar mass fraction as the systems most favourable to disk onset under these assumptions.

We evaluate this diagnostic against our statistical sample in Figure 13, colouring each snapshot of a galaxy by the $v/\sigma$ value for star-forming gas, stars younger than 10 Myr, or the full stellar population. The snapshots form a continuous distribution in the $R_{1/2}/R_{\rm v}$–$M_{\star}/M_{\rm v}$ plane, with no clear structural bimodality. Young stars show a clear kinematic gradient across this plane: rotational support is more common in galaxies with a larger stellar mass fraction and a more centrally concentrated stellar component. The full stellar population shows a similar but weaker trend, consistent with the older stars retaining a broader range of orbital motions.

Image not extracted

In contrast, the star-forming gas exhibits a distinct kinematic behavior. Numerous systems with rotation-supported young stars host dispersion-dominated star-forming gas—most notably at low stellar masses where the galaxy-to-halo size ratio, $R_{1/2}/R_{\rm v}$, remains low. Conversely, rotation-supported gas is predominantly found in higher-mass systems that have developed a centrally concentrated stellar component. This difference matters because young stars inherit the resolved velocities of their parent gas but do not retain its thermal or unresolved turbulent pressure support in our $v/\sigma$ measurement. A coherent orbiting clump or interaction-triggered burst can therefore produce rotation-supported young stars even when the star-forming gas as a whole is not a settled disk. Young-star rotation alone may consequently overstate disk formation in low-mass galaxies. The full stellar population provides a complementary, longer-term measure, but its response is diluted by pre-existing stars and their accumulated orbital history.

This raises the separate question of long-term disk survival. Dekel et al. (2020) argue that disks become long-lived only above $M_{\rm v} \simeq2 \times10^{11}\,M_{\odot}$, where the merger timescale exceeds the orbital time ($t_{\rm mer} > t_{\rm orb}$). At our median stellar-to-halo mass ratio of $M_{\star}/M_{\rm v} \simeq0.02$, this halo mass corresponds to $M_{\star} \simeq4 \times10^{9}\,M_{\odot}$, close to the transition mass identified empirically in § 5. This coincidence suggests that the characteristic stellar mass for disk onset in our sample may also approximate the stellar-mass threshold for long-term disk survival, although agreement in mass alone does not establish a common cause.

The diagnostic of Eq. (6) closely resembles the analytical boundary derived by Semenov (2026) (purple line in their Fig. 5). Using a more complete treatment, their framework places the turnover on the critical line $v_{\rm c,tot}(R_{1/2,\rm b}) = v_{\rm v}$, whereas our implementation differs in two practical respects. First, we use the stellar rather than the total circular speed, which makes our diagnostic more conservative. Second, we evaluate the rotation curve at the stellar half-mass radius instead of the baryonic half-mass radius. We use the stellar half-mass radius because it can be measured from stellar-mass maps and is less sensitive to the adopted galaxy aperture in our analysis (Stevens et al. 2014). Stellar dominance within $R_{1/2}$ makes the stellar and total circular speeds similar there, but it does not make the stellar and baryonic half-mass radii identical, since an extended gas component can substantially increase the latter.

The approximation $v/\sigma\sim v_{\rm c}/v_{\rm v}$ fails when $v_{\rm rot} \ll v_{\rm c}$ or $\sigma\ll v_{\rm v}$ (Eq. (5)). An observational example for the first case is found in the stellar dynamical models of CALIFA galaxies by Kalinova et al. (2017), where galaxies with sharply peaked circular-speed curves are slowly rotating and pressure-supported systems (their Fig. 9). A central peak in the inferred circular speed therefore need not correspond to strong stellar rotation. Gas can reach $v_{\rm rot} \simeq v_{\rm c}$ more readily than a pre-existing collisionless stellar population because shocks and radiative cooling dissipate non-circular motions, provided that the gas has coherent angular momentum. However, collisionality alone does not ensure this outcome, since pressure support, continuing inflow, and feedback can prevent the gas from settling. A deep central potential thus facilitates rotation without guaranteeing it.

The second case, $\sigma\ll v_{\mathrm{v}}$, usually applies at later epochs, after a dynamically cold, thin disk has settled. Übler et al. (2019) measure a decline in the mean intrinsic ionized-gas dispersion from $\sim45$ km/s at $z \sim2.3$ to $\sim30$ km/s at $z \sim0.9$, well below the virial velocity of the host haloes, and Tamburro et al. (2009) find a characteristic H i dispersion of $10 \pm2$ km/s at the optical radius of 11 nearby disk galaxies. These measurements trace different gas phases and should not be equated directly, but together they illustrate a dynamically colder regime in which $\sigma\ll v_{\mathrm{v}}$ and rotational support can be high even without a centrally peaked rotation curve. At early epochs, before disk formation, gas velocity dispersion can instead be comparable to the virial velocity, as shown in simulations (Goldner et al. 2026). In our sample $\sigma_{\mathrm{SFgas}} \simeq0.6 v_{\mathrm{v}}$ both before and after the transition (panel (b) of Figure 12).

Because disk formation criteria are sensitive to the chosen definition of $v/\sigma$ and underlying theoretical assumptions, simulations diverge on the precise criteria. Specifically, young stellar tracers systematically exhibit stronger rotational support than the underlying star-forming gas, and standard kinematic approximations such as $v_{\mathrm{rot}} \simeq v_{\mathrm{c}}$ and $\sigma\simeq v_{\mathrm{v}}$ could break down. In the ARTEMIS zoom-in simulations of Milky Way-mass halos, Dillamore et al. (2024) find that stellar disks generally form at a halo virial mass threshold of $\simeq(6 \pm3) \times10^{11}\,M_{\odot}$, while the circular velocity profile is still rising between one and three stellar half-mass radii (their Figs. 4 and 5). Similarly, using the FIRE-2 suite, Myrtaj et al. (2026) show that their Romeo system develops a young stellar disk while the circular-velocity slope diagnostic remains positive. They attribute this disk formation to the cessation of baryonic mass sloshing relative to the center of mass within the host potential, noting that the velocity slope becomes negative only several gigayears after spin-up (their Fig. 20). In MAGE, we instead identify disk formation using star-forming gas weighted by SFR. Under this metric, reaching a sufficiently high stellar mass alongside an established, centrally concentrated stellar core provides a more robust criterion for disk formation. We note that the disks formed here are kinematically thick rather than thin Milky Way analogs, as we discuss in § 6.5, meaning that $\sigma$ remains comparable to the halo virial velocity.

Impact of feedback strength on galaxy sizes and kinematics

We now assess the sensitivity of our results to the adopted feedback model, with particular emphasis on photoionization heating. Fig. 14 compares four simulations of the MAGE07 halo at $z = 9$ that differ only in the strength of early stellar feedback from photoionization: highphot, fiducial (our MAGE07 run), lowphot, and nophot, in order of decreasing photoionization strength. The model parameters are described in Andalman et al. (2025). For each run, we compare the galaxy-to-halo size ratio ($R_{1/2}/R_{\mathrm{v}}$), the stellar-to-halo mass ratio ($M_{\star}/M_{\mathrm{v}}$), and the rotational support of stars younger than 10 Myr.

The figure shows that varying photoionization strength produces systematic shifts in the normalized size–mass plane. From highphot to nophot, the stellar-to-halo mass ratio, $M_{\star}/M_{\mathrm{v}}$, increases from $\simeq0.02$ to $\simeq0.08$, while the galaxy-to-halo size ratio, $R_{1/2}/R_{\mathrm{v}}$, decreases from $\simeq0.024$ to $\simeq0.009$. The two models with the weakest photoionization, lowphot and nophot, have similar stellar-to-halo mass ratios and small normalized sizes, with size ratios of approximately 0.010 and 0.009, respectively.

Image not extracted

The rotational support of the newly formed stars also increases as photoionization weakens: $v/\sigma$ rises monotonically from $\simeq1.1$ in highphot to $\simeq1.7$ in fiducial, $\simeq2.3$ in lowphot, and $\simeq2.8$ in nophot. Rotational support continues to increase from lowphot to nophot, although their mass and size ratios change comparatively little. The stellar surface-density maps illustrate the corresponding structural differences. With strong photoionization, the stellar body is extended and diffuse, spanning most of the 2 kpc field, whereas with weak or no photoionization the stellar distribution is concentrated in a compact, high-surface-density core surrounded by a low-density envelope and satellite clumps. The compact cores in the models with weak or no photoionization also appear elongated in projection, consistent with their higher $v/\sigma$.

The joint trend in size, mass, and kinematics is consistent with the evolutionary sequence identified in § 5. Weaker photoionization allows gas to cool and collapse more efficiently in the centre, which enhances central star formation, deepens the central potential, and raises $v_{\mathrm{c,tot}}(R_{1/2})$ relative to $v_{\mathrm{v}}$. The weaker-photoionization runs evolve further along the compaction sequence of Figure 12, towards the compact, more rotationally supported regime. Stronger photoionization, acting together with supernova feedback, may favour more spatially extended star formation and larger stellar half-mass radii. While highphot remains extended and @page 20/30 starts_mid=1 ends_mid=1 a pre-existing collisionless stellar population because shocks and radiative cooling dissipate non-circular motions, provided that the gas has coherent angular momentum. However, collisionality alone does not ensure this outcome, since pressure support, continuing inflow, and feedback can prevent the gas from settling. A deep central potential thus facilitates rotation without guaranteeing it.

The second case, $\sigma\ll v_{\mathrm{v}}$, usually applies at later epochs, after a dynamically cold, thin disk has settled. Übler et al. (2019) measure a decline in the mean intrinsic ionized-gas dispersion from $\sim45$ km/s at $z \sim2.3$ to $\sim30$ km/s at $z \sim0.9$, well below the virial velocity of the host haloes, and Tamburro et al. (2009) find a characteristic H i dispersion of $10 \pm2$ km/s at the optical radius of 11 nearby disk galaxies. These measurements trace different gas phases and should not be equated directly, but together they illustrate a dynamically colder regime in which $\sigma\ll v_{\mathrm{v}}$ and rotational support can be high even without a centrally peaked rotation curve. At early epochs, before disk formation, gas velocity dispersion can instead be comparable to the virial velocity, as shown in simulations (Goldner et al. 2026). In our sample $\sigma_{\mathrm{SFgas}} \simeq0.6 v_{\mathrm{v}}$ both before and after the transition (panel (b) of Figure 12).

Because disk formation criteria are sensitive to the chosen definition of $v/\sigma$ and underlying theoretical assumptions, simulations diverge on the precise criteria. Specifically, young stellar tracers systematically exhibit stronger rotational support than the underlying star-forming gas, and standard kinematic approximations such as $v_{\mathrm{rot}} \simeq v_{\mathrm{c}}$ and $\sigma\simeq v_{\mathrm{v}}$ could break down. In the ARTEMIS zoom-in simulations of Milky Way-mass halos, Dillamore et al. (2024) find that stellar disks generally form at a halo virial mass threshold of $\simeq(6 \pm3) \times10^{11}\,M_{\odot}$, while the circular velocity profile is still rising between one and three stellar half-mass radii (their Figs. 4 and 5). Similarly, using the FIRE-2 suite, Myrtaj et al. (2026) show that their Romeo system develops a young stellar disk while the circular-velocity slope diagnostic remains positive. They attribute this disk formation to the cessation of baryonic mass sloshing relative to the center of mass within the host potential, noting that the velocity slope becomes negative only several gigayears after spin-up (their Fig. 20). In MAGE, we instead identify disk formation using star-forming gas weighted by SFR. Under this metric, reaching a sufficiently high stellar mass alongside an established, centrally concentrated stellar core provides a more robust criterion for disk formation. We note that the disks formed here are kinematically thick rather than thin Milky Way analogs, as we discuss in § 6.5, meaning that $\sigma$ remains comparable to the halo virial velocity.

Impact of feedback strength on galaxy sizes and kinematics

We now assess the sensitivity of our results to the adopted feedback model, with particular emphasis on photoionization heating. Fig. 14 compares four simulations of the MAGE07 halo at $z = 9$ that differ only in the strength of early stellar feedback from photoionization: highphot, fiducial (our MAGE07 run), lowphot, and nophot, in order of decreasing photoionization strength. The model parameters are described in Andalman et al. (2025). For each run, we compare the galaxy-to-halo size ratio ($R_{1/2}/R_{\mathrm{v}}$), the stellar-to-halo mass ratio ($M_{\star}/M_{\mathrm{v}}$), and the rotational support of stars younger than 10 Myr.

The figure shows that varying photoionization strength produces systematic shifts in the normalized size–mass plane. From highphot to nophot, the stellar-to-halo mass ratio, $M_{\star}/M_{\mathrm{v}}$, increases from $\simeq0.02$ to $\simeq0.08$, while the galaxy-to-halo size ratio, $R_{1/2}/R_{\mathrm{v}}$, decreases from $\simeq0.024$ to $\simeq0.009$. The two models with the weakest photoionization, lowphot and nophot, have similar stellar-to-halo mass ratios and small normalized sizes, with size ratios of approximately 0.010 and 0.009, respectively.

Image not extracted

The rotational support of the newly formed stars also increases as photoionization weakens: $v/\sigma$ rises monotonically from $\simeq1.1$ in highphot to $\simeq1.7$ in fiducial, $\simeq2.3$ in lowphot, and $\simeq2.8$ in nophot. Rotational support continues to increase from lowphot to nophot, although their mass and size ratios change comparatively little. The stellar surface-density maps illustrate the corresponding structural differences. With strong photoionization, the stellar body is extended and diffuse, spanning most of the 2 kpc field, whereas with weak or no photoionization the stellar distribution is concentrated in a compact, high-surface-density core surrounded by a low-density envelope and satellite clumps. The compact cores in the models with weak or no photoionization also appear elongated in projection, consistent with their higher $v/\sigma$.

The joint trend in size, mass, and kinematics is consistent with the evolutionary sequence identified in § 5. Weaker photoionization allows gas to cool and collapse more efficiently in the centre, which enhances central star formation, deepens the central potential, and raises $v_{\mathrm{c,tot}}(R_{1/2})$ relative to $v_{\mathrm{v}}$. The weaker-photoionization runs evolve further along the compaction sequence of Figure 12, towards the compact, more rotationally supported regime. Stronger photoionization, acting together with supernova feedback, may favour more spatially extended star formation and larger stellar half-mass radii. While highphot remains extended and only marginally rotationally supported, both lowphot and nophot host compact galaxies with well-developed young stellar disks and $M_\star\simeq6 \times10^{9}\,M_\odot$, above the disk-formation mass of $\simeq5 \times10^{9}\,M_\odot$. According to the criterion in Eq. (6), both ratios change in ways that favour disk formation as photoionization weakens: at fixed halo mass, the stellar-to-halo mass ratio increases while $R_{1/2}/R_{\rm v}$ decreases.

These experiments vary photoionization while keeping the SN prescription fixed, so they do not establish the relative importance of photoionization and SN feedback. They demonstrate that physical galaxy sizes, stellar masses, and young-star rotational support depend on the feedback prescription. This sensitivity is not specific to cosmic dawn. Zoom-in simulations of Milky Way-mass haloes at lower redshift have shown that the strength and implementation of stellar feedback change galaxy sizes at fixed halo mass by factors of a few and determine whether an extended disk forms at all (Sales et al. 2010; Übler et al. 2014; Agertz & Kravtsov 2016). Galaxy sizes have accordingly served as a calibration target for subgrid feedback in large-volume simulations (Crain et al. 2015).

Comparison to previous theoretical work

Compact sizes at the low-mass end

At low stellar masses, the difference between our sizes and those of other simulations is not a modest offset in normalization but a change in the nature of the objects. At $M_\star\simeq10^{8}\,M_\odot$, our galaxies have half-mass radii of $\simeq20$–$30$ pc, resembling dense star clusters rather than the extended, faint stellar bodies formed at the same mass in THESAN-ZOOM, FirstLight, and Lumina (Figure 4 and Figure 5). The FIRE-2 galaxies of Ma et al. [85] are similarly extended: at $z = 10$, their projected half-mass radii are $\simeq0.3$–$0.5$ kpc at $M_\star\simeq10^{7}$–$10^{8}\,M_\odot$, approximately half those in THESAN-ZOOM but an order of magnitude larger than ours.² Other suites, however, do form compact systems. In SERRA (Pallottini et al. 2022), central galaxies at $z = 7.7$ with $M_\star\lesssim10^{9}\,M_\odot$ also have small reported sizes, with dust-attenuated UV half-light radii $\lesssim100$ pc; these cannot be compared directly with our intrinsic sizes. The recent VESTRAL simulations also produce compact galaxies at high redshift (Nelson et al., in preparation; Ravishankar et al., in preparation), and future comparisons will allow us to assess similarities and differences between their results and ours.

Such differences are commonly attributed to the sub-grid star formation and feedback models. These models certainly affect galaxy sizes: varying effective photoionization feedback changes $R_{1/2}/R_{\rm v}$ of MAGE07 by a factor of $\simeq3$ (Section 6.3). It is nevertheless difficult to attribute an offset of more than an order of magnitude to this effect alone. SERRA provides a direct test because it also uses RAMSES, but adopts an otherwise independent set of sub-grid models, including a Schmidt–Kennicutt law applied to molecular hydrogen with a fixed efficiency of 10 per cent per free-fall time; SN, stellar-wind, and radiation-pressure feedback with blast-wave energy dissipation; non-equilibrium chemistry with KROME; and on-the-fly radiative transfer with RAMSES-RT. Despite these differences, its central galaxies at $z = 7.7$ have sizes within a factor of a few ours.

Numerical factors beyond the sub-grid models may therefore matter. The AGORA code comparison at $z \geq10$ (Kim et al. 2026) finds that RAMSES produces noticeably higher stellar masses and star-formation rates at $z \simeq13$–$15$ than both the particle-based codes and ENZO, regardless of halo mass, before converging with the other codes to within $\sim1$ dex by $z \sim10$, once the haloes collapse and the star-forming gas becomes dense in all codes. They attribute the early boost to the refinement strategy of RAMSES, which extends the maximally refined region around dense cells to avoid abrupt resolution jumps between dense and neighbouring lower-density regions. The resulting expansion of the highly refined star-forming region allows RAMSES reach $n_{\rm H} \sim3 \times10^{3}\,{\rm cm}^{-3}$ in every halo to have a efficient local star cluster formation (Grudić et al. 2018; Menon et al. 2024), whereas the other codes reach such densities only in the most massive systems. This suggests that numerical factors, in addition to the sub-grid prescriptions, may contribute to higher star-forming gas densities. Such effects could also occur in MAGE and SERRA, depending on their refinement settings. Denser star-forming gas favours more concentrated star formation and a deeper central potential, promoting the formation of star clusters. However, whether the refinement strategy also affects stellar sizes remains to be tested, as Kim et al. (2026) do not report these measurements.

Sample selection effects may also play a role. Because both MAGE and SERRA target massive galaxies at cosmic dawn via cosmological zoom-in simulations, our galaxy size measurements are likely influenced by the overdense environments sampled by these zoom-in simulations, particularly at the low-mass end. According to the FFB scenario (Dekel et al. 2023), galaxies exceeding the threshold denoted by the red line in Figure 1 are expected to be more compact and host bound star clusters. At $z > 12$, this threshold corresponds to a halo mass of $M_{\rm v} \simeq10^{10}\,M_\odot$ and a stellar mass of a few times $10^{8}\,M_\odot$, which could naturally account for the observed compactness of low-mass systems at these early epochs. However, our low-mass galaxies at $z = 9$–$11$ remain comparably compact, and the three target galaxies lying below the FFB threshold at $z = 9$ show no qualitative structural differences from those above it, suggesting that the transition between the FFB and standard regimes is gradual rather than sharp. In addition, the large-size tail in SERRA is driven predominantly by non-central galaxies undergoing mergers or having recently merged. Because our sample is restricted to central galaxies and explicitly excludes systems with multiple dense stellar cores, obvious major mergers are filtered out (§ 2.2). We note, however, that including merger systems primarily inflates the size scatter without substantially altering the median relation.

Turnover of the size–mass relation

The turnover of the intrinsic stellar size–mass relation occurs at different masses across simulation suites. In THESAN-1, Shen et al. (2024) find a maximum in the stellar half-mass radius near $M_\star\sim10^{8}\,M_\odot$ at $z \gtrsim6$. In FirstLight, Cataldi et al. (2026) identify a negative slope during a compaction phase spanning $M_\star\sim10^{8.5}$–$10^{9.5}\,M_\odot$ over $5 \leq z \leq9$. Negative intrinsic slopes have also been reported at lower redshifts for $M_\star\gtrsim10^{9}\,M_\odot$ in BLUETIDES (Marshall et al. 2022), ILLUSTRIS-TNG (Genel et al. 2018), and FLARES (Roper et al. 2022). Our turnover near $M_\star\simeq5 \times10^{9}\,M_\odot$ therefore occurs at a higher mass, with an additional difference in size normalization (Figure 4). We note that these comparisons concern intrinsic stellar half-mass radii. A turnover in the intrinsic half-mass size–mass relation need not persist in the relation measured using dust-attenuated UV half-light radii (§ 6.6.1).

The proposed mechanisms all link central gas inflow to stellar mass buildup. Cataldi et al. (2026) attribute the negative slope to compaction triggered by inflows that strengthen as the galactic centre becomes baryon-dominated, although they do not identify the physical mechanism that triggers these inflows. Shen et al. (2024) instead associate rapid compaction in massive galaxies with internal disk instability. In THESAN-ZOOM, McClymont et al. (2025) find repeated central starbursts followed by inside-out quenching and renewed size growth, while the population retains a positive intrinsic size–mass relation. Thus, central bursts and individual contraction episodes do not necessarily produce a turnover in the population relation. Their cumulative effect on the central mass distribution, and the subsequent spatial distribution of star formation, determine whether size growth resumes or remains suppressed.

Our feedback experiments demonstrate that changes in the size normalization cannot be interpreted independently of changes in stellar mass. Increasing photoionization feedback in the same halo produces a less massive, more extended stellar body at $z = 9$ (Section 6.3), moving the galaxy towards the larger sizes found in the comparison samples. However, this also reduces the stellar contribution to the central potential, since $v_{\mathrm{c,star}}^2(R_{1/2}) = GM_{\star}/(2R_{1/2})$, which can delay the development of a central potential deep enough to overcome the effects of feedback. A shift towards lower stellar mass and larger radius induced by stronger feedback therefore does not imply that the turnover itself moves to lower mass. The four feedback variants establish the sensitivity of one halo’s structure to photoionization, but do not measure how the population’s transition mass changes with feedback strength.

A further consideration is that MAGE does not model black-hole growth or AGN feedback. At lower redshift, zoom-in simulations that switch AGN feedback on and off in the same haloes find that it enlarges massive galaxies and removes the flattening or turnover of the size–mass relation at high mass. In RAMSES zoom-ins of massive haloes evolved to $z = 0$, AGN jets turn overly massive, compact, and rotation-dominated central galaxies into extended, dispersion-dominated systems that follow the observed size–mass relation (Dubois et al. (2013)). Choi et al. (2018) find that, without AGN feedback, the 30 massive galaxies have present-day sizes approximately one-fifth of those observed at fixed stellar mass and follow a shallow growth track in the size–mass plane, whereas mechanical and radiative AGN feedback quenches in-situ star formation, so that subsequent gas-poor mergers and the expulsion of central gas puff up the stellar body along the observed relation. In FIRE, Cochrane et al. (2023) re-simulate a $M_{\star} \sim10^{11}\,\mathrm{M}_{\odot}$ galaxy at $z \simeq2.3$ and show that, without AGN winds, a rapid central starburst shrinks the half-mass radius, whereas mechanical winds lower the central gas surface density, suppress the compact starburst, and hold the half-mass radius approximately constant. In FIRE-3 runs evolved to $z = 0$, the most massive galaxies are an order of magnitude too compact without AGN feedback, and the stronger feedback model recovers the observed size–mass relation (Byrne et al. (2024)). Even at $z \simeq6$, bursty AGN feedback in simulated quasar hosts expands compact bulges by an order of magnitude by $z \simeq3$ through fluctuations of the gravitational potential (van der Vlugt & Costa [161]). The processes identified in these studies affect the physical ingredients underlying our turnover: efficient central star formation, the retention of gas in a deepening central potential, and the resulting rapid central stellar buildup.

At cosmic dawn, however, a direct test in one of our haloes finds that AGN feedback has a comparatively weak effect. Sunseri et al. (2026) introduced a central black hole into the MAGE07 halo, systematically varying the seed mass ($10^{4}$–$10^{5}\,\mathrm{M}_{\odot}$), the accretion limit (one to three times the Eddington rate (as well as unlimited accretion)), the thermal AGN feedback coupling efficiency, and the strength of stellar feedback. Because stellar feedback maintains a turbulent, multiphase nuclear ISM, the central black hole remains starved for roughly half the time after the onset of star formation. Even when AGN feedback becomes self-regulating, the star formation history, the final stellar mass at $z \simeq9$, and the thick, partly rotationally supported disk remain essentially unchanged relative to the black-hole-free baseline. Furthermore, while THESAN-1 includes AGN feedback and still exhibits a negative slope, THESAN-ZOOM lacks AGN feedback yet displays a positive relation. AGN feedback alone therefore does not guarantee a positive slope. If early black holes grow more rapidly than assumed in these numerical experiments and cosmological simulations, as suggested by the proposed AGN interpretation of GN-z11 (Maiolino et al. (2024)), stronger feedback could potentially flatten or erase the negative high-mass slope. Consequently, the turnover mass reported here should be interpreted strictly within a stellar-feedback-only framework, although the parameter explorations of Sunseri et al. (2026) and relevant simulations suggest that the impact of omitting AGN feedback may be modest for haloes assembling as rapidly as ours at $z \gtrsim9$.

Disk formation

Comparing disk formation across simulations requires distinguishing the emergence of coherent rotation from subsequent settling into a dynamically cold, thin disk. In the FIRE-2 simulations, Yu et al. [154] and McCluskey et al. (2024) identify successive phases in which stars form first in a dispersion-dominated protogalaxy, then in a turbulent rotating disk, and finally in a colder disk. The transition to thin-disk formation is associated with steadier star formation and virialization of the inner CGM (Stern et al. 2021; Hafen et al. 2022; Gurvich et al. 2023). The inner CGM virializes when the cooling time of shocked gas at $\sim0.1R_{\mathrm{v}}$ exceeds the local free-fall time, so that cool, supersonic inflows and outflows with large pressure fluctuations give way to a hot, subsonic, quasi-static phase near the virial temperature. Accreting gas can then align its angular momentum with that of the galaxy as before joining the galaxy. The gas velocity dispersion drops sharply and can fall well below the virial velocity.

The connection between central stellar mass buildup and rotational support in MAGE is qualitatively consistent with the controlled experiments of Hopkins et al. (2023). They show that a centrally concentrated mass profile promotes disk formation by providing a well-defined dynamical centre and suppressing large-scale motions that disrupt coherent rotation. In TNG50, Semenov et al. (2024) find that disk spin-up accompanies bulge growth, the inflow of prealigned gas, outflow recycling, and hot-halo formation. Their analysis does not establish which process causes the transition. Here, “bulge” refers to a central mass concentration and does not necessarily imply dispersion support: Luo et al. (2026) find that, at $z = 3$, the central structures formed through compaction in massive TNG50 galaxies are mostly rotation-supported nuclear disks and bars, which become dispersion dominated only at $z \lesssim1$. The analytical model of Semenov (2026) links disk formation to the steepening of the potential as baryons accumulate and predicts a corresponding change in size regulation. However, as discussed in § 6.2, the ARTEMIS and FIRE-2 results show that a centrally peaked rotation curve need not precede the formation of a stellar disk (Dillamore et al. 2024; Myrtaj et al. 2026).

Simulations resolving the cold interstellar medium also demonstrate that disks can form well before the disk-formation epochs reported for Milky Way progenitors in many earlier studies. Semenov et al. (2025b) find that a disk forms at $z \sim6$–$7$ in an early-forming Milky Way analogue. Disk formation is aided by rapid central mass growth and reduced gas pressure support. Companion experiments by the same authors show that a turbulence-dependent star-formation efficiency can allow both bursty star formation at early times and disk survival (Semenov et al. 2025a). At $z > 7$, Mayer et al. (2025) find compact, self-gravitating gas disks that fragment into massive star-forming clumps in an overdense environment. These studies provide precedents for early rotating systems.

An alternative route to the early formation of massive disks invokes cosmology rather than internal galaxy physics. Shen et al. (2026b) present cosmological hydrodynamic simulations with Early Dark Energy, in which the modified pre-recombination expansion history accelerates early structure formation and raises the number density of stellar and gaseous disks by $\simeq0.5$ dex at $z \simeq6$–$7$. This enhancement primarily reflects the higher abundance of massive galaxies. At $z \lesssim3$, the model predictions approach those of $\Lambda$CDM$.

Other simulations using different subgrid models for star formation and stellar feedback also suggest an association between compaction and disk formation. In THESAN-1, Shen et al. (2024) find that more massive, compact systems have a larger fraction of their gas mass in the disk component, consistent with our results. By contrast, in THESAN-ZOOM, McClymont et al. (2025) report a positive intrinsic size–mass relation, while Wang et al. (2025) find that galaxies in the same suite have not yet developed well-defined disks. Thus, even if our size turnover and disk formation are not linked by a direct causal relationship, the broader association between compactness and disk formation appears in several simulation frameworks.

Observational implications

The size evolution and the connection between central stellar mass buildup and disk formation in MAGE can be tested through the joint structural, star-formation, and kinematic properties of high-redshift galaxy observations.

Size–mass and size–magnitude relations

M24 find a positive UV size–mass relation, $R_{\mathrm{e}} \propto M^{0.19 \pm0.03}$, and relatively weak size evolution, $R_{\mathrm{e}} \propto(1 + z)^{-0.4 \pm0.2}$, across $5 < z < 14$. Their sample also contains a compact subpopulation ($\simeq13$ per cent of their galaxies) with intense star formation, which is qualitatively consistent with our positive low-mass branch and weak evolution over $9 \lesssim z \lesssim15$ (Figure 4 and Figure 5). The COSMOS-Web sample of Yang et al. (2025), which extends to lower redshift, instead yields stronger evolution of optical sizes over $2 < z < 10$, with $R_{\mathrm{e}} \propto(1 + z)^{-1.21 \pm0.05}$ at a reference mass of $5 \times10^{10}\,M_{\odot}$. The size–mass slopes are poorly constrained above $z \sim8$. They also report very compact blue sources that remain unresolved at $z > 5$, even after excluding X-ray AGN and little red dots. At later epochs, Chen et al. (2026b) find that the rest-frame optical size–mass relation of star-forming galaxies at $0.5 < z < 6$ flattens above a nearly redshift-independent pivot mass of $M_{\star} \sim10^{10}\,M_{\odot}$, within a factor of two of our turnover mass of $\simeq5 \times10^{9}\,M_{\odot}$. They interpret this as evidence that galaxy size growth tracks halo growth below the pivot, whereas central mass buildup decouples the size growth of more massive galaxies from halo growth, in line with the saturation of stellar birth radii and the decline of $R_{1/2}/R_{\mathrm{v}}$ in MAGE (Figure 4 and Figure 10 and § 6.1).

In the size–magnitude relation shown in Figure 6, all eight C iv-strong sources are labelled, and all but GHZ4 have half-light radii $\lesssim100$–$150$ pc, with most sources remaining unresolved (Roberts-Borsani et al. (2026)). They lie on or just above the upper envelope of our simulated distribution, and because most of their reported sizes are upper limits, their intrinsic sizes may be smaller still. Their strong rest-frame UV lines, including C iv and N iv] in GHZ2 and MoM-z14 (Castelló-Mor et al. (2024); Naidu et al. (2026)), are indicative of hard ionizing radiation fields from intense, compact star formation and/or active galactic nucleus (AGN) activity (Harikane et al. (2025); Roberts-Borsani et al. (2026)). The evidence for a stellar origin of these lines is discussed in § 6.7.2. The C iv-weak sources, in contrast, range from similarly compact, unresolved systems to extended galaxies with $R_{\mathrm{e}} \simeq1$ kpc. Compactness is therefore not always accompanied by strong lines. U37126 is a clear example: it is similarly compact ($R_{\mathrm{e}} \simeq61$ pc; see footnote 1) and extremely blue ($\beta\simeq-2.9$) but shows no detectable nebular emission. Marques-Chaves et al. (2026) attribute this absence to a near-unity Lyman-continuum escape fraction after feedback cleared gas from the stellar core of U37126. Interestingly, GN-z11, the brightest source in the sample, is classified as C iv-weak in the prism data, but its spectrum shows strong N iv], N iii], and C iii] emission and a possible Cygni C iv feature (Bunker et al. (2023)). Its extreme compactness and strong emission lines may reflect an intense, weakly obscured central starburst (Chen et al. 2026a; Nakane & Ouchi 2026) and/or a luminous unresolved AGN that enhances the central UV surface brightness (Maiolino et al. (2024); Scholtz et al. (2024)).

Although our simulation results match the observed compact systems, we do not reproduce the more extended population indicated by the gray bands in Figure 4 and Figure 6. Dust attenuation represents one of the primary systematic uncertainties in comparing intrinsic UV sizes with observed size relations. Including dust attenuation in THESAN-1 (Shen et al. 2024) and FIRSTLIGHT (Ceverino et al. (2026)) changes their intrinsically negative size–mass slopes to flat or positive slopes, bringing them into closer agreement with observations. This behavior arises because intrinsically compact, massive galaxies frequently host concentrated central dust columns. Preferential attenuation of the luminous core consequently expands the emergent UV half-light radius. By contrast, THESAN-ZOOM shows a negligible effect of dust on UV half-light radii, a result mirrored by the ASTRID simulation for rest-frame optical sizes at lower redshifts ($z = 3$–$6$; LaChance et al. (2025)). The extent to which dust alters UV sizes therefore depends on the spatial distribution of dust and its survival under stellar feedback. Rapid dust destruction in star-forming environments minimizes attenuation (Garaldi et al. (2026)).

Observationally, the compact systems with strong emission lines identified by Roberts-Borsani et al. (2026) have sizes overlapping those in our simulations and blue UV continuum slopes, $\beta< -2.5$. These blue slopes suggest little UV attenuation, potentially because of radiation-pressure-driven dust clearing or spatial segregation of dust from star-forming complexes (Ferrara et al. (2023); Ziparo et al. (2023)). In contrast, typical systems with weak emission lines exhibit redder slopes between $\beta\simeq-2.5$ and $-1.5$, leaving room for moderate size expansion. For MAGE, preliminary radiative transfer calculations suggest that dust attenuation can increase the effective UV half-light radius by a factor of a few in our massive galaxies (Nakazato et al., in preparation). In the size–magnitude plane, central obscuration increases the half-light radius while simultaneously diminishing the total UV luminosity. Given that Kim et al. (2026) report significant UV attenuation even at a dust-to-metal ratio of 0.01, dust shifts galaxies along both axes of Figure 6. However, the mean UV continuum slope of the $z > 9$ population becomes bluer with increasing redshift, changing from $\langle\beta\rangle\simeq-2.2$ at $z \simeq9.5$ to $\langle\beta\rangle\simeq-2.6$ at $z \simeq11.5$, approaching the theoretical limit for dust-free stellar populations (Topping et al. (2024); Cullen et al. (2024); Morales et al. (2024); Austin et al. (2025)). In this interpretation, dust-induced size inflation may be more important for a dust-rich minority at the lower-redshift end of our sample, whereas observed UV half-light radii at $z > 11$ should be close to their intrinsic values. Future work will determine the net shift in the MAGE size–magnitude relation and test whether such attenuation can be reconciled with the blue UV slopes at $z \gtrsim11$. Moreover, at low stellar masses, where the discrepancy between our results and observations is more severe, the effect of dust is limited. Another factor is bias arising from limited observational resolution. The NIRCam F200W PSF has a FWHM of $\simeq0.066^{\prime\prime}$, corresponding to $\simeq250$–$300$ pc at $z = 9$–$12$, whereas the UV half-light radii of our galaxies with $M_\star\lesssim10^{8.5}\,\mathrm{M}_\odot$ are $R_{\rm e} \simeq10$–$40$ pc, an order of magnitude smaller. Profile fitting with a PSF model can recover sizes below the FWHM, but sizes far below it are poorly constrained ((Ono et al. 2023; M24)), so such galaxies enter observed size–mass and size–magnitude relations with size upper limits or measurements scattered upwards towards the resolution floor. The compact subpopulation of M24 is defined partly by this limit, as it includes gravitationally lensed sources with $R_{\rm e}$ below half the PSF FWHM. Detecting additional compact systems at the low-mass end could steepen the inferred size–mass relation and bring it into closer agreement with our results.

Although higher-resolution observations may yield smaller median sizes, our simulations produce no galaxies with intrinsic UV half-light radii approaching $\simeq1$ kpc, even among the most massive systems. On the observational side, a fraction of seemingly extended sources may be close merging pairs that cannot be cleanly deblended in projection. On the simulation side, our sample selection preferentially targets overdense regions as mentioned in § 6.5.1. Because our low-mass systems are close neighbors of massive zoom-in halos rather than typical field galaxies, they assemble earlier and reside in halos that are more concentrated at fixed mass than those probed by blank-field surveys. Consequently, our suite underrepresents the broader field population, which could artificially compress the predicted scatter in galaxy size. A rigorous quantitative comparison will require forward-modeled mock observations with matched prescriptions for dust attenuation, PSF convolution, surface brightness noise, source selection, and profile fitting, which we defer to future work.

In-situ star formation and merger-driven growth

Recent JWST studies infer different contributions from mergers to early stellar mass growth. For a mass-selected JADES sample at $z \simeq3$–$9$, Puskás et al. (2025) estimate that stellar accretion through major mergers contributes approximately 3–13 per cent of the total growth, with in-situ star formation dominating. Using close pairs at $4.5 < z < 11.5$, Duan et al. (2025) infer a larger contribution: major-merger stellar accretion amounts to $71 \pm25$ per cent of the growth through star formation. Both analyses infer accretion statistically from pair counts and merger time-scales, rather than directly identifying the birthplaces of observed stars. Puskás et al. (2025) note that their merger rates agree with those of Duan et al. (2025), and suggest that different estimates of the specific SFR help explain the discrepant growth fractions.

Our stellar birth-site analysis provides a complementary constraint on where the accreted and locally formed populations reside. Within $R_{1/2}$, the median in-situ fraction in each mass group is $\gtrsim0.7$. The four high-mass systems span $f_{\rm in\mathchar-situ} \simeq 0.8$–$1$, while ex-situ stars contribute more strongly in the outskirts ({xref}Fig. 8). Moreover, the total half-mass radius closely follows that of the in-situ component ({xref}Fig. 9), supporting a dominant contribution from the in-situ population without isolating newly formed stars from subsequent orbital redistribution. Such high in-situ fractions are consistent with THESAN-ZOOM, where {citet}Shen et al. (2026a)` find that the fraction of stellar mass formed within the virial radius of the main progenitor is close to unity for halos with $M_{\rm vir} \lesssim10^{10}\,\mathrm{M}_\odot$ at $z \simeq3$–$10$ and remains above $\simeq0.6$ even in more massive haloes at lower redshift, where merger-delivered stars become noticeable. This supports the qualitative picture of predominantly in-situ central growth.

However, in-situ dominance does not imply that mergers are dynamically negligible. Gas-rich interactions can directly funnel gas or induce torque-driven inflows that trigger rapid starbursts within the central galaxy on timescales of a few tens of megayears, as documented in both observations ((Harikane et al. 2025; Kiyota et al. 2025)) and cosmological simulations ((Nakazato et al. 2024)). Because stars formed during these bursts are born from the central gas reservoir, they are classified as in situ despite their merger-driven origin. The MAGE04 system exemplifies this behavior, maintaining a high central in-situ fraction even after undergoing recent major mergers (Figure 7).

Early disk formation

Observational evidence for early disks spans structural and kinematic diagnostics. Ono et al. (2025) find that galaxies at $z \sim10$–$16$ have a log-normal size distribution, stellar-to-halo size ratios of $R_{\rm e}/R_{\rm v} \simeq0.015$, a mean Sérsic index close to unity, and a nearly uniform axis-ratio distribution, which they interpret as consistent with disk formation. At $4 < z < 7.7$, Mirzoyan et al. (2026) use a denoising neural-network morphological classifier to classify $\simeq34$ per cent of 100 JWST/NIRCam galaxies as disk-like across a range of stellar masses. At $z = 7.43$, Baker et al. (2025) identify a compact stellar core with a half-light radius of $\sim80$ pc surrounded by a star-forming component with a half-light radius $\simeq400$ pc. The specific SFR increases with radius, providing evidence for inside-out growth around an already dense centre. These observations support the possibility that compact central buildup and subsequent disk growth can occur early, as in our simulations.

Structural similarities alone do not establish the degree of rotational support. Gas kinematics reveal a range of dynamical states at these epochs. Parlanti et al. (2023) find that ALMA [C ii] and [O iii] observations of 22 star-forming galaxies at $z > 5$ are consistent with turbulent rotating disks, although the available resolution does not uniquely distinguish disks from closely separated mergers. Using JWST/NIRSpec spectroscopy, de Graaff et al. (2024) infer both rotation- and dispersion-dominated configurations among six low-mass galaxies at $5.5 < z < 7.4$. A broader population study by Danhaiye et al. (2025) finds that $34 \pm25$ per cent of H$\alpha$-selected galaxies with $\log M_\star/M_\odot< 10^{10}$ have $v/\sigma_0 > 1$ in both the $z \sim4.5$ and $z \sim5.5$ bins, where $\sigma_0$ is the intrinsic ionized-gas velocity dispersion. In the Twister-Z5 galaxy at $z = 5.4$, Nelson et al. (2024) infer rapid ionized-gas rotation, with $v = 305 \pm70$ km s$^{-1}$ and $v/\sigma\sim2$–$3$. The increase in H$\alpha$ equivalent width with radius is consistent with inside-out growth. Together, these observations show that ordered rotation can coexist with substantial gas turbulence.

Dynamically cold gas disks are also observed, primarily through [C ii] observations, which trace dynamically colder gas than H$\alpha$ observations ((Kohandel et al. 2024)). Rizzo et al. (2021) infer $v/\sigma\sim7$–$15$ from ALMA [C ii] kinematics of strongly lensed dusty star-forming galaxies at $z \sim4$–$5. Lelli et al. (2026) find regular rotation in all 16 massive star-forming galaxies ($M_\star\simeq10^{10.5}$–$10^{12}\,\mathrm{M}_\odot$) of the TRICEPS sample at $z = 4$–$5$, observed in [C ii] at $\simeq0.8$ kpc resolution. Regular rotation is present even though half of the galaxies belong to pairs, suggesting that rotational support was widespread among massive systems by that epoch. Mock [C ii] observations of massive TNG50 galaxies at the same redshifts yield $v/\sigma\simeq2$–$3$ in most cases ((He et al. 2026)), below the observed values, so reproducing such cold disks remains a challenge for current simulations. At still earlier epochs, Rowland et al. (2024) establish the presence of a dynamically cold disk in REBELS-25 at $z = 7.31$ using ALMA [C ii] observations, while Rowland et al. (2026) find that [C ii] and [O iii] 88 µm trace a common rotating structure in REBELS-25, despite the galaxy’s clumpy rest-frame UV morphology. Similarly, the strongly lensed Cosmic Grapes galaxy at $z = 6.072$ contains at least 15 star-forming clumps embedded in a rotating gas disk (Fujimoto et al. 2025). Such systems show that clumpy emission can coexist with ordered rotation and does not by itself imply a merger.

Tentative kinematic evidence extends into the redshift range probed by MAGE. For GN-z11 at $z = 10.6$, Xu et al. [151] model the C iii] velocity gradient as disk rotation and infer $v/\sigma= 2.83^{+1.82}_{-1.41}$, while noting that an outflow can also explain the observations. At $z = 14.18$, Scholtz et al. (2025) report a tentative detection of a velocity gradient at $3\sigma$ significance in [O iii] 88 µm emission from JADES-GS-z14-0. Their disk model gives $v/\sigma> 2.5$. These candidates suggest that appreciable rotational support may emerge very early, but higher-resolution observations are needed to test whether the observed velocity gradients arise from disk rotation.

A proposed evolutionary sequence and its observational implications

Figure 15 combines the trends of § 4 and § 5 and the interpretation developed above into a proposed evolutionary sequence, ordered by stellar mass, together with its possible observational counterparts.

Physical sequence

The sequence is organized by the relative magnitudes of three velocity scales: the halo virial speed $v_{\mathrm{v}}$, the total circular speed at the half-mass radius $v_{\mathrm{c,tot}}(R_{1/2})$, and the characteristic SN-feedback scale $v_{\mathrm{SN}} \simeq100\ \mathrm{km\,s^{-1}}$ (§ 5). We interpret the trends as four approximate regimes, correspond to the gas-escape, halo-recycling, galactic-recycling, and young-disk regimes in Figure 15. Panel labels below refer to Figure 12.

(i) Gas escape and compact, cluster-like stellar bodies, $M_{\star} \lesssim10^{8}\ M_{\odot}$: $v_{\mathrm{v}} < v_{\mathrm{SN}}$, so SN-driven winds can remove a substantial fraction of the gas from the halo itself (Dekel & Silk 1986), and dark matter dominates the enclosed mass at $R_{1/2}$ (panel (a)). The star-forming gas and the young stars remain dispersion-dominated, with median $v/\sigma\simeq0.5$–$1$ (panel (d)) and only transient excursions to higher values (Figure 11). Stars form predominantly outside the half-mass radius, at $R_{\mathrm{birth}} \simeq2R_{1/2}$, with birth radii varying by more than an order of magnitude between galaxies, consistent with the chaotic assembly of the low-mass progenitors (Figure 10). The galaxies remain compact, with half-mass radii only a few times the spatial resolution of $\simeq10\ \mathrm{pc}$.

(ii) Halo-scale retention and outward star formation, $10^{8}\ M_{\odot} \lesssim M_{\star} < 5 \times10^{8}\ M_{\odot}$: $v_{\mathrm{v}} > v_{\mathrm{SN}}$ while $v_{\mathrm{c,tot}}(R_{1/2}) < v_{\mathrm{SN}}$, and dark matter still dominates the enclosed mass at $R_{1/2}$ (panel (a)). SN-driven outflows are expected to be retained by the halo but not by the galaxy. In this interpretation, ejected gas may mix with high-angular-momentum CGM gas before returning to the galactic outskirts. The measured median birth radius is $R_{\mathrm{birth}} \simeq2R_{1/2} \simeq0.02\ R_{\mathrm{v}}$ (Figure 10), and $R_{1/2}$ grows steeply with stellar mass (panel (c)). Fresh inflow appears to dominate the central gas supply: near the lower-mass boundary, the measured outflow rate across $R_{1/2}$ amounts to $\lesssim0.1$ of the estimated net inflow rate across $R_{\mathrm{gal}}$.

(iii) Galaxy-scale retention and central baryonic buildup, $5 \times10^{8}\ M_{\odot} \lesssim M_{\star} < 5 \times10^{9}\ M_{\odot}$: $v_{\mathrm{c,tot}}(R_{1/2}) > v_{\mathrm{SN}}$, and stars dominate the enclosed mass at $R_{1/2}$. We interpret this transition as a shift from halo-scale to galactic-scale recycling. The outflow rate across $R_{\mathrm{gal}}$ falls below that across $R_{1/2}$ (panel (e)). This ordering is consistent with failed outflows and recycling within the galaxy that could sustain the buildup of central stellar mass. Size growth slows, and the central gas fraction declines as the stellar contribution grows.

(iv) Size turnover, young-disk onset, and gas-poor centres, $M_{\star} \gtrsim5 \times10^{9}\ M_{\odot}$: $v_{\mathrm{c,tot}}(R_{1/2}) > v_{\mathrm{v}}$, driven by the stellar component, while $v_{\mathrm{rot,SFgas}}/v_{\mathrm{c,tot}}(R_{1/2})$ rises from $\simeq0.4$ to $\simeq0.7$ and $\sigma_{\mathrm{SFgas}}$ continues to decrease with mass (panel (b)). The rapid central stellar buildup at $M_{\star} \simeq5 \times10^{9}\ M_{\odot}$ coincides with the turnover of $R_{1/2}$ (panel (c)) and a rise in rotational support: both the star-forming gas and stars younger than 10 Myr reach $v/\sigma\gtrsim1.5$, while the full stellar population remains dispersion-dominated (panel (d)). The outflow rate across $R_{1/2}$ becomes comparable to the estimated net inflow rate across $R_{\mathrm{gal}}$, suggesting that recycled gas may contribute as much as fresh accretion to the central supply (panel (e)). The central star-formation fraction remains at $\simeq0.4$ up to the turnover. At the highest masses sampled, $M_{\star} \gtrsim2.5 \times10^{10}\ M_{\odot}$, MAGE01 shows a decline in this fraction, while gas contributes only a few per cent of the central enclosed mass (panels (a) and (e)). This decline, the low central gas fraction, and the younger outskirts (Figure 8) suggest the onset of inside-out quenching.

Observational counterparts

The sequence offers a possible connection between the compact and extended populations introduced in § 1. Compact galaxies with strong high-ionization lines may be observed near central burst peaks, preferentially during the galactic-recycling phase, when the dust cycle is high and youngest stars dominate the UV light from a region of a few tens of parsecs (Harikane et al. 2025; Roberts-Borsani et al. 2026). For the $z > 10$ sample of Roberts-Borsani et al. (2026), individual and composite spectra support burst star formation as the primary explanation for the spectral diversity. The strong C IV emitters have the bluest UV slopes ($\beta< -2.5$) and the highest inferred star-formation surface densities, consistent with young, concentrated bursts.

GN-z11 at $z = 10.6$ provides a plausible example of this phase and may be at the transition to disk formation. Its UV light is dominated by an unresolved component embedded in more extended emission, with a combined half-light radius of $R_{\mathrm{e}} = 64 \pm20\ \mathrm{pc}$ (Tacchella et al. 2023). Stellar-wind features indicate a recent burst with an age of $\lesssim3\ \mathrm{Myr}$ (Chen et al. 2026a). Stellar models better reproduce its UV absorption and emission diagnostics than the AGN models considered (Nakane & Ouchi 2026), and narrow rest-frame optical lines support vigorous star formation (Álvarez-Márquez et al. 2025). An accreting black hole has nevertheless been proposed based on the reported [Ne IV] and C II* detections, high gas densities, and a fast ionized outflow (Maiolino et al. 2024; Scholtz et al. 2024). The tentative rotation inferred by Xu et al. [151] suggests that central star formation may coexist with rotational support. Gas recycling may contribute to the metal enrichment invoked in models of its chemical abundances (Citro et al. 2026). Its estimated stellar mass of order $10^{9}\ M_{\odot}$ (Tacchella et al. 2023) falls naturally within this galactic recycling regime, reinforcing our proposed framework.

Compact galaxies with weak or undetected lines may be the same systems caught between bursts, when SN feedback has cleared the gas and nebular lines fade within $\sim10\ \mathrm{Myr}$ while the UV continuum uum remains concentrated in the stellar core. These lulls should be longer below $5 \times10^{8}\,M_{\odot}$, where expelled gas returns through the halo, than above it, where galactic recycling restores the supply more quickly (§ 6.1). Weak lines can also occur soon after a burst if its gas has already been removed: U37126 has a very young, extremely blue stellar population and a half-light radius of $\simeq61$ pc but shows no detectable nebular lines. This absence is attributed to a near-unity Lyman-continuum escape fraction after gas clearing ((Marques-Chaves et al. 2026)). If longer replenishment times produce a lower burst duty cycle, lower-mass galaxies may more frequently show weak or undetected nebular lines. Quantifying this prediction requires detailed emission-line modelling and the observational selection function.

More extended galaxies may correspond to later disk growth around a less active core, after inside-out quenching suppresses central star formation (§ 6.1 and § 6.6.1). Their light profiles are consistent with this possibility: typical star-forming galaxies at $z \simeq1$–$6$ have median Sérsic indices of $n \simeq1$–$1.5$ (Shibuya et al. 2015), while the luminous galaxies of Ono et al. (2025) at $z \sim10$–$16$ have $n = 1.3 \pm0.6$ and typical sizes of 100–300 pc. The size and axis-ratio distributions reported by Ono et al. provide additional morphological evidence for disks (§ 6.6.3). However, this mapping is an extrapolation: our runs end at $z \simeq9$, shortly after the turnover, and the simulated half-mass radii do not re-grow within the sampled range. Renewed size growth, seen after compaction in FIRSTLIGHT (Cataldi et al. 2026) and VELA (Lapiner et al. 2023), is therefore shown in Figure 15 as an expectation rather than a result of our simulations.

Mergers offer another route to extended morphologies. In FOR-EVER22, the stellar component becomes temporarily larger during a major merger and shrinks after coalescence and angular-momentum loss (Yajima et al. 2022; Ono et al. 2025) likewise invoke mergers or tidal interactions for some of their largest galaxies. Our main-branch tracks show transient enlargement during the major merger of MAGE04 (Figure 3), but the statistical sample excludes snapshots with multiple spatially separated stellar components (§ 2.2). Such interacting systems therefore have no counterpart in our size–mass relation.

The associations apply most directly when star formation powers the compact emission. The observed spectral diagnostics do not uniquely establish the source of nebular ionization: for GHZ2 at $z = 12.34$, UV line ratios are consistent with both AGN photoionization and low-metallicity star formation, although the limit on [Ne IV]/N IV] favours stellar photoionization in some model grids (Castellano et al. 2024). More generally, UV line ratios from metal-poor star formation at high ionization parameter overlap with those of AGN, and narrow lines do not exclude an obscured nucleus. Because our simulations omit black-hole growth (§ 6.8), an unresolved AGN contribution must be separated before comparing observed sizes with our intrinsic stellar UV sizes. Such a contribution can enhance the central surface brightness and reduce the measured half-light radius (§ 3.3).

The proposed counterparts suggest distinct observational tests. Later disks should preferentially be massive systems with ordered rotation and older central populations. Compact weak-line galaxies should retain young or intermediate-age populations in gas-cleared cores, while interacting systems should show multiple components or tidal features. Resolved stellar ages, emission-line maps, and kinematics spanning the compact and extended populations would test whether they include different phases of the proposed sequence. In GN-z11, such measurements could distinguish a brief central burst from sustained stellar buildup. These associations are tentative and need not identify a unique evolutionary phase. The limitations of the sample and physical model are discussed in § 6.8.

Caveats

Sample and selection. The main limitations are the small number of independent haloes and the selection of massive haloes in overdense regions. The low-mass galaxies in these zoom regions inhabit the overdense environments around the massive target haloes and do not constitute a representative field sample. This environmental difference may contribute to their compactness relative to observed low-mass galaxies. To measure the structure of individual stellar systems, we impose a cut on the outer stellar-mass fraction to exclude active merger phases, during which sizes are temporarily inflated. Consequently, this selection introduces a bias against systems with multiple spatially separated stellar components. The high-mass turnover is driven primarily by four systems, and only MAGE01 contributes at the largest stellar masses. These systems cross the transition mass only shortly before $z = 9$, so the phase after disk onset is sampled briefly and the half-mass radii do not re-grow within the simulated range. The identification of extended observed systems with later disk growth in § 6.7 is therefore an extrapolation. Repeated snapshots are correlated, and different progenitor tracks can share an environment or a descendant, so regression errors that treat snapshots as independent may underestimate the full sampling uncertainty, particularly for the high-mass slope. A larger sample is needed to determine the intrinsic scatter and any redshift dependence of the characteristic stellar mass inferred from the combined structural and kinematic trends.

Feedback model. Absolute sizes and stellar masses depend on the adopted feedback model. Four feedback variants of one halo demonstrate sensitivity to photoionization but do not establish how population slopes or transition masses depend on feedback strength. The simulations also omit black-hole growth and AGN feedback, as well as radiation pressure and stellar winds. These processes could modify the central stellar buildup and gas retention identified here, particularly in the most massive systems. The omission of AGN feedback bears most directly on the high-mass turnover: at lower redshift, simulations run with and without AGN feedback show that it can suppress compact central starbursts, enlarge massive galaxies, and remove the flattening or turnover of the size–mass relation at high mass ([161] (Dubois et al. 2013; Choi et al. 2018; Cochrane et al. 2023; Byrne et al. 2024)). The negative high-mass slope and the turnover mass should therefore be regarded as predictions of a model with stellar feedback alone (§ 6.5.2).

Radiation pressure. The neglect of radiation pressure deserves particular emphasis. Our early feedback model is limited to photoionization heating (§ 2), yet analytic and radiation-hydrodynamic studies of massive star-cluster formation find that radiation pressure, rather than thermal pressure of photoionized gas, becomes the dominant pre-SN feedback channel at the high gas surface densities and cluster masses characteristic of starburst environments (Krumholz & Matzner 2009; Murray et al. 2010; Fall et al. 2010; Kim et al. 2018; Krumholz et al. 2019). Radiation-hydrodynamic simulations of cluster-forming clouds under the dense ($\Sigma\gtrsim10^{3}\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-2}$), metal-poor conditions expected at cosmic dawn find that direct radiation pressure on dust, rather than reprocessed infrared radiation, limits the star-formation efficiency and launches outflows from such clouds (Menon et al. 2022, 2024, 2025), with Lyman-$\alpha$ radiation pressure potentially contributing where dust is scarce (Ferrara et al. 2025; Menon & Smith 2026). The omission of radiation pressure may therefore cause MAGE to underestimate the strength of early stellar feedback in the dense, compact star-forming regions that host the rapid central star formation identified here. By analogy with the photoionization variants of Section 6.3, stronger early feedback could lower stellar masses and enlarge radii at fixed halo mass. However, the same cloud-scale studies find that the star-formation efficiency remains high at low dust-to-gas ratios. The net effect on the transition mass identified here is therefore uncertain.

Numerical resolution. The smallest stellar radii approach the numerical resolution, making the sizes of the compact low-mass galaxies of § 6.5.1 potentially sensitive to resolution. Resolution tests are needed to assess the sensitivity of the main results.

Correspondence to observed quantities. The quantities presented here are intrinsic properties of the simulated galaxies and are not directly equivalent to observational measurements. The stellar half-mass radii are not directly observable, and the UV sizes we present here omit dust attenuation and instrumental effects such as PSF convolution and noise. Likewise, neither the star-forming gas disks nor the stellar disks in § 5 correspond directly to observationally inferred early disks. Such observed systems are typically traced by cool or ionized gas via emission lines such as [C II], [O III], H$\alpha$, and C III]. Measurements of these systems depend on aperture choices, beam smearing, and the adopted velocity-dispersion definitions. Consequently, our gas or stellar $v/\sigma$ values should not be directly equated with observationally inferred gas $v/\sigma$. A more direct comparison with observations therefore requires additional post-processing: radiative-transfer calculations for the dust-attenuated UV light, emission-line modelling of the gas kinematics, and mock observations with matched apertures, PSF, source selection, and fitting procedures.

Gas flows and stellar redistribution. Shell fluxes do not distinguish returning gas from fresh accretion and do not account for the motion of the aperture boundary. They therefore support a recycling interpretation without measuring individual recycling times or the complete mass balance of a moving region. The birth-radius statistics identify where stars form but do not quantify how subsequent orbital redistribution changes galaxy sizes. Tracking gas trajectories would constrain recycling, while following fixed stellar cohorts would quantify stellar redistribution.

Summary

We have investigated galaxy size evolution and disk formation over $9 \lesssim z \lesssim15$ using MAGE, a new suite of 11 cosmological zoom-in simulations performed with RAMSES at $\sim10\,\mathrm{pc}$ spatial resolution using the star-formation and feedback model of Andalman et al. (2025). The target galaxies span $M_{\star} \simeq2 \times10^{8}$–$3 \times10^{10}\,M_{\odot}$ at $z = 9$. Our main conclusions are as follows.

(i) Galaxy sizes increase steeply with stellar mass at low masses and evolve weakly with redshift at fixed mass. Near $M_{\star} \simeq5 \times10^{9}\,M_{\odot}$, the stellar half-mass size–mass relation turns over, while the UV half-light size–mass relation becomes approximately flat (Figure 4). The ratio $R_{1/2}/R_{\mathrm{e}}$ reaches $\simeq0.02$ near the transition and subsequently declines. Our intrinsic UV sizes overlap the compact subpopulation of M24 and the compact end of the observed size–magnitude distribution (Figure 4 and Figure 6), but do not reproduce the more extended observed galaxies. (ii) In-situ stars dominate the central stellar mass, and the half-mass radius of the in-situ component closely tracks that of the full stellar population. At $z = 9$, the in-situ fractions within $R_{1/2}$ are $\gtrsim0.7$ in all but MAGE08, which experienced a major merger within the preceding 50 Myr (Figure 7 and Figure 8). Characteristic stellar birth radii remain near $0.02 R_{\rm v}$ below $M_\star\simeq5 \times10^9\,M_\odot$ and saturate around 400 pc above this mass, suggesting a transition from size growth that tracks halo growth to stronger internal regulation (Figure 9 and Figure 10).

(iii) Comparing the characteristic SN-feedback scale $v_{\rm SN} \simeq100\,\mathrm{km\,s^{-1}}$ (Dekel & Silk 1986) with the halo and central velocity scales identifies two characteristic masses below the size turnover: the virial velocity exceeds $v_{\rm SN}$ near $M_\star\simeq10^8\,M_\odot$, and the total circular speed at $R_{1/2}$ exceeds it near $M_\star\simeq5 \times10^8\,M_\odot$. The outflow rates measured across $R_{1/2}$, $R_{\rm gal}$, and $R_{\rm v}$ show a radial ordering consistent with increasing confinement of outward-moving gas as the potential deepens. The SFR within $R_{1/2}$ remains at $\simeq0.3$–$0.55$ of the galaxy-wide SFR until $M_\star\gtrsim2.5 \times10^{10}\,M_\odot$, where this fraction declines in MAGE01, accompanied by a decline in $v_{\rm c,gas}$. Together with the younger outskirts, these changes suggest the onset of inside-out quenching (Figure 8 and Figure 12).

(iv) The rapid central stellar buildup drives $v_{\rm c,tot}(R_{1/2})$ above $v_{\rm v}$ near $M_\star\simeq5 \times10^9\,M_\odot$ and coincides with the size turnover. Rotational support rises rapidly in the star-forming gas and young stars, with both the gas and stars younger than 10 Myr reaching $v/\sigma\simeq1.5$–$3$ at $z \simeq9$ (Figure 11 and Figure 12). For the star-forming gas, the rise in $v/\sigma$ is driven primarily by increasing rotation speed: $v_{\rm rot,SFgas}/v_{\rm c,tot}(R_{1/2})$ increases from $\simeq0.4$ to $\simeq0.7$, while its effective dispersion stays near $0.6v_{\rm v}$. The full stellar population remains dispersion-dominated. Together with the preceding central buildup, this joint structural and kinematic transition shares key features with wet compaction events at cosmic noon.

(v) Varying photoionization feedback in the same halo changes both size and stellar mass by factors of a few. Weaker photoionization feedback produces a more massive, compact stellar body and increases young-star rotational support from $v/\sigma\simeq1.1$ to $\simeq2.8$, indicating a correlation between central concentration and rotation in these experiments (Fig. 14).

Figure 15 and § 6.7 summarize the proposed evolutionary sequence, ordered by stellar mass, in which the deepening central potential increasingly confines SN-driven outflows. At cosmic dawn, GN-z11 combines evidence for compact central star formation (Chen et al. 2026a; Nakane & Ouchi 2026) with tentative evidence for rotational support from its C iii] velocity gradient (Xu et al. [151]). This makes it a promising candidate for exploring the connection between central stellar buildup and early disk formation. The inferred transition masses depend on the adopted feedback model and limited sample, while direct observational tests require matched modelling of dust attenuation and kinematic tracers, which we leave to future work. The Nancy Grace Roman Space Telescope’s tiered High-Latitude Wide-Area Survey could identify tens to thousands of galaxies at $z \sim13$ (Yung et al. [155]), depending on the bright end of the UV luminosity function (Harikane et al. 2022). Such galaxies could provide counterparts to our massive systems, and the brightest could be targets for JWST and ALMA follow-up of their sizes, star-formation profiles, and kinematics.

Acknowledgements

We dedicate this paper to the memory of Avishai Dekel, our co-author, mentor and friend. May his memory be a blessing. We thank Fangzhou Jiang, Hyeonyong Kim, Sharon Lapiner, Jinning Liang, Rongrong Liu, Yurina Nakazato, Dylan Nelson, Anirudh Ravis-hankar, Xuejian Shen, and Jonathan Stern for helpful discussions. We thank the VESTRAL and LUMINA teams for sharing results ahead of publication, and the FIRSTLIGHT team for providing dust-free UV half-light radii. This material is based upon work supported by the National Science Foundation (NSF) and the U.S.-Israel Binational Science Foundation (BSF) under Award Number 2406558 and Award Title “The Origin of the Excess of Bright Galaxies at Cosmic Dawn”. ZY and NM acknowledge support from BSF grant 2022736 (from the NSF-BSF program) and from BSF grant 2022281. The authors gratefully acknowledge the Scientific Computing Core at the Flatiron Institute for providing the high-performance computing resources and storage that enabled the research results reported within this paper. The Flatiron Institute is supported by the Simons Foundation. Artificial intelligence tools assisted with code development, language editing, literature verification, and the rendering of the artwork in the conceptual schematic of Figure 15, whose content, layout, and annotations were specified and verified by the authors. The authors retain responsibility for the scientific analysis, interpretation, conclusions, and accuracy of the manuscript.

Data Availability

The simulation data underlying this article will be shared on reasonable request to the corresponding author.

This paper has been typeset from a TeX/LaTeX file prepared by the author.