© The Authors 2026

1. Introduction

Sizes and internal kinematics constitute fundamental proxies for understanding the processes involved in the formation of galaxies, and those shaping their structures. These proxies are the result of many internal and external multi-scale physical mechanisms, such as galaxy mergers, instabilities, gas accretion and outflows, feedback processes, and star formation. However, the precise connection between the size and kinematics of the stellar component and these mechanisms is still poorly known (e.g. Hopkins et al. 2023, and references therein).

The Sloan Digital Sky Survey has provided a large statistical sample of galaxies in the local Universe, measuring reliable data on their properties (e.g. Kauffmann et al. 2003). These data, together with those from the Hubble Space Telescope (HST), have allowed to investigate fundamental relations that help to unveil the aforementioned connection, such as that between stellar size and mass (Shen et al. 2003; Simard et al. 2011; Cappellari 2013; van der Wel et al. 2014; Mowla et al. 2019a,b). This relation presents a positive slope at low redshifts (z ≲ 2), indicating that galaxies grow inside-out.

The picture is not so clear at higher redshifts. The HST has revealed a population of massive galaxies ($M^{★}\sim {10}^{10}{\mathrm{M}}_{\odot }$) at z ∼ 7 that are more compact than their counterparts at lower redshifts (e.g. Bruce et al. 2012; van der Wel et al. 2014; Allen et al. 2017; Bouwens et al. 2017; Yang et al. 2021). The James Webb Space Telescope (JWST) has unveiled an unexpected population of red galaxies at z ∼ 7 − 9, with stellar masses of M^★ ≳ 10¹⁰ M_⊙, extremely small effective radii (⟨r_e⟩≈150 pc), high star formation rate (SFR) surface densities (Σ_SFR > 10 M_⊙ yr⁻¹ kpc⁻²), and positive evidence of inside-out growth (Baggen et al. 2023; Baker et al. 2025; Morishita et al. 2024; Ormerod et al. 2024; Ward et al. 2024; Miller et al. 2025; Yang et al. 2025). Also Tacchella (2025) has found compact galaxies with sizes of the order of 100 pc in the JADES survey. Interestingly, Baggen et al. (2023) report an inverted (i.e. negative slope) stellar size–mass relation, whereas Morishita et al. (2024) argue that, while the stellar size–UV luminosity relation is inverted, the size–mass relation is normal. The question of how massive galaxies become compact and the inverted size–mass relation (if confirmed) arises should therefore be addressed.

Theoretical models suggest that the normal size–mass relation arises because the gas infalling to the central regions conserves its angular momentum, forming a disc that grows in an inside-out way (Kravtsov et al. 1997; Mo et al. 1998). Star formation proceeds therefore at larger radii as the mass of the system increases. High-resolution cosmological hydrodynamical simulations confirm this picture (Genel et al. 2014; Pedrosa & Tissera 2015; Pillepich et al. 2018; Popping et al. 2022), and indeed are capable of reproducing the relation at low redshifts (Furlong et al. 2016).

Challenges appear in galaxies at high redshifts. Numerical studies of galaxies at the Cosmic Dawn, with different baryonic implementation and numerical resolution (BLUETIDES, FLARES IX, THESAN; Ni et al. 2020; Lovell et al. 2021; Kannan et al. 2022), obtain an intrinsic size–mass relation with a maximum at a stellar mass of between 10⁸ to 10⁹ M_⊙. At lower masses galaxies are influenced by the strength of feedback-driven outflows, resulting in growing sizes as mass increases. In contrast, massive galaxies are affected by disc instability that triggers compaction (i.e. a decrease in the size of the stellar distribution¹), inverting the relation. Also Roper et al. (2023) find a transition between extended and irregular objects to compact ones via efficient centralised cooling that produces a high SFR in their cores. However, they find a normal size–mass relation and argue that the observed inversion of the size–luminosity trend is due to the effect of dust.

Several works explore the mechanisms behind size changes, either theoretically (Dekel & Burkert 2014) or using cosmological simulations (Zolotov et al. 2015; Lapiner et al. 2023). They present different versions of the ‘wet compaction’ scenario, in which the dissipative collapse of high-redshift, gas-rich, turbulent discs produces compact star-forming systems (‘blue nuggets’) at z ≲ 4 − 5. These discs are fed by cold streams or mergers, and undergo instabilities that drive gas inwards. Compaction triggers central gas depletion and quenching, making galaxies evolve into compact spheroids (‘red nuggets’) by z ∼ 2. Quenching proceeds inside-out, and thus galaxies resume expansion as star formation moves outwards into a ring (Dekel et al. 2020). At low redshift, blue nuggets are rare due to reduced gas content. According to Lapiner et al. (2023), wet compaction occurs near a critical stellar mass of ∼10¹⁰ M_⊙. Driven by angular momentum loss, this transition reshapes their structure, composition, and kinematics – from diffuse to compact, dark matter to baryon central dominance, and supernova (SN)- to active galactic nucleus (AGN)-driven feedback. Also, Shen et al. (2024) find that low-mass galaxies (M^★ ∼ 10^7 − 9 M_⊙) at z > 3 do not undergo a compaction–re-expansion process, but rapid size fluctuations produced by the competition between feedback-driven gas outflows and cold inflows. These fluctuations preserve the normal size–mass relation, enlarging its dispersion instead of changing its slope.

Given the dissimilar results about the high-redshift size–mass and size–luminosity relations found by both observations and various simulations, a thorough investigation is necessary to elucidate the involved processes. However, numerical resolution of a few hundred parsecs from large-scale cosmological simulations such as ILLUSTRISTNG (Pillepich et al. 2018), THESAN (Shen et al. 2024), FLARES (Roper et al. 2022), BLUETIDES (Marshall et al. 2022), and SIMBA (Wu et al. 2020) makes hard to resolve internal structures at high redshift.

High-resolution simulations such as THESAN-ZOOM (McClymont et al. 2025) or FIREBOX (Feldmann et al. 2023) are the present workhorses used to overcome this limitation. In this paper, we use the FIRSTLIGHT (FL) simulation suite (Ceverino et al. 2017, 2018, 2019, 2021, 2024), which achieves a high spatial resolution of ∼10 pc and provides a large sample size for a robust analysis. We used snapshots of 169 simulated galaxies, each at 45 redshifts in the range 5.25 − 9 (a total of 7605 snapshots). We aim to study the growth of the stellar distribution of galaxies at Cosmic Dawn, and the physical mechanisms responsible for the high-mass galaxy compaction and the inversion of the mass–size relation.

The structure of this paper is as follows. Sect. 2 provides a brief synopsis of the FL simulations used throughout this work. Sect. 3 describes our results on the evolution of the sizes of galaxies, whereas Sect. 4 compares them to those obtained by other works, discusses their relevance to the understanding of galaxy evolution in the early Universe, and presents our conclusions.

2. The simulation

This paper uses a sub-sample (SS) from the cosmological zoom-in simulation suite described in detail in Ceverino et al. (2017). The FL project has been capable of predicting the mechanisms driving galaxy morphology (Nakazato et al. 2024; Ceverino et al. 2021). They have also produced a galaxy population that matches a wide range of observations at high redshifts such as star-formation histories (Ceverino et al. 2018, 2024), rest-frame UV/optical absolute magnitudes and colours, mass-metallicity (Langan et al. 2020) and mass-dust relations (Mushtaq et al. 2023), and optical emission lines (Ceverino et al. 2019, 2021; Nakazato et al. 2023).

The simulations were performed with the ART code (Kravtsov et al. 1997; Kravtsov 2003; Ceverino & Klypin 2009), which accurately follows the evolution of a gravitating N-body system and the gas dynamics using an Eulerian adaptive mesh refinement (AMR) approach. Besides gravity and hydrodynamics, the code incorporates many astrophysical processes relevant to galaxy formation as subgrid physics that include gas cooling due to atomic hydrogen and helium, metal and molecular hydrogen cooling, photoionization heating by a constant cosmological UV background with partial self-shielding (Haardt & Madau 1996), star formation, and feedback (thermal, kinetic, and radiative).

Stellar particles are created in gas cells at numerical densities above a threshold of 1 cm⁻³ and temperatures below 10⁴ K. The code implements a stochastic star formation model that yields the empirical Kennicutt–Schmidt law (Schmidt 1959; Kennicutt 1998; Ceverino & Klypin 2009). The simulations track the metals released from SNe Ia and SNe II, using yields from Woosley & Weaver (1995).

In addition to thermal energy feedback the simulations use radiative feedback, modelled as a non-thermal pressure added to the total gas pressure in regions where ionising photons from massive stars are produced and trapped. It uses moderate trapping of infrared photons (details in Ceverino et al. 2014). The code also includes the injection of momentum from the (unresolved) expansion of gaseous shells from SNe and stellar winds (Ostriker & Shetty 2011). More details can be found in Ceverino et al. (2017), Ceverino & Klypin (2009), and Ceverino et al. (2010, 2014). The effect of the different feedback models is discussed in Ceverino et al. (2014, 2023). AGN feedback is not included in these simulations. It may have an effect on quenching star formation at very high stellar masses, M^★ > 10¹⁰ M_⊙, where black holes are massive enough to influence their host galaxies (Nelson et al. 2019; Lapiner et al. 2023).

The galaxy sample consists of 169 galaxies tracked in the redshift interval 5.25 < z < 9, and whose haloes have a maximum circular velocity range between 50 and 300 km s⁻¹. This covers a halo mass range between 10⁹ and a few times 10¹¹ M_⊙, and more than four orders of magnitude in stellar mass, 10⁶ < M^★/M_⊙ < 10^10.5. We stored a total of 45 snapshots for each galaxy, with a time spacing of 7 − 10 Myr².

The target haloes were initially selected using low-resolution N-body-only simulations of three cosmological boxes with sizes 10, 20, and 40 h⁻¹ Mpc, assuming WMAP5 cosmology with matter density, baryon density, Hubble, and density fluctuation parameters Ω_m = 0.27, Ω_b = 0.045, h = 0.7, and σ₈ = 0.82, respectively (Komatsu et al. 2009). Initial conditions for the selected haloes were then generated using a standard zoom-in technique (Klypin et al. 2011), with a much higher resolution. The resulting resolution of dark matter particles is m_DM = 10⁴ M_⊙, with a minimum stellar particle mass of 100 M_⊙. The maximum spatial resolution is always between 8.7 and 17 proper parsecs (a comoving resolution of 109 pc after z = 11).

3. Results

To represent the size of FL simulated galaxies, we used the half-mass radius of the distribution of stellar particles, $r_{\mathrm{hm}}^{★}$. This is defined as the radius of a sphere that encloses half of the stellar mass found within 0.15 r₂₀₀, with r₂₀₀ the usual definition of the virial radius (i.e. the radius of a spherical volume with a mean overdensity of 200). We adopted this quantity because it is less influenced by geometric effects than 2D-projected sizes, and also by tracer biases such as those affecting the sizes computed from the light in specific bands. For all galaxies we find at least 1000 stellar particles within $r_{\mathrm{hm}}^{★}$ at any redshift, which makes our analysis numerically robust and allows us to achieve a well-sampled stellar mass distribution.

In Fig. 1 we show the distribution of $r_{\mathrm{hm}}^{★}$. Its peak radius and width evolve to larger values with decreasing redshift, displaying the hierarchical growth of galaxies as they are assembled. To subtract the bulk of this effect from the general size evolution, we also present the distribution of the stellar-to-halo size ratio ${\xi }^{★}=r_{\mathrm{hm}}^{★}/r_{200}$. A bimodality is already apparent; most galaxies belong to the peak at log ξ^★ ≈ −1 (hereafter regular galaxies). This subpopulation shows a very weak evolution with redshift, in the sense that haloes grow faster than the distribution of stars, which is expected because the latter is regulated by feedback whereas the former are not. The second peak at log ξ^★ ≲ −1.6 indicates the existence of a subpopulation of compact galaxies. This peak grows as redshift decreases, indicating a continued production of these objects. The fact that our sample comprises only those galaxies that can be followed through the whole redshift interval implies that the compact subpopulation grows at the expense of the regular one. In other words, there are galaxies that become more compact as time elapses. The median value of ξ^★ over the full redshift range (5.25 < z < 9) is ⟨ξ^★⟩ = 0.07, similar to the value of 0.08 reported by Ma et al. (2018) for FIRE zoom-in simulated galaxies at comparable epochs (6 < z < 10). However, Ma et al. (2018) do not find evidence of bimodality in their sample. This is likely due to the fact that the majority of their galaxies are low-mass systems with M^★ < 10⁸ M_⊙. Furthermore, they report an overall expansion in their galaxy population, in contrast with our findings.

Fig. 1. Distribution of the half-mass stellar radius (left panel) and the stellar-to-halo size ratio (right panel), for 5.25 < z < 9 (colour coded). The dashed vertical line indicates the median value of ξ★ over the full sample and redshift range (⟨ξ★⟩ = 0.07). A growing subpopulation of compact galaxies is seen; the dotted vertical line at log ξ★ = −1.6 roughly illustrates the separation.

Fig. 2 plots the relation between the stellar mass, M^★, of galaxies and the virial mass, M₂₀₀, of their haloes. We define the former as the stellar mass within a radius enclosing 83% of the baryonic mass of the galaxy (called the optical radius) and the latter as the total mass within r₂₀₀. FL galaxies follow a tight relation in the M₂₀₀ − M^★ plane, with a scatter of ∼1 dex, assumed to be driven by the different mass accretion histories (Moster et al. 2018). On the one hand (top panel) this relation is consistent with the results of semi-empirical models of Tacchella et al. (2018) and Behroozi et al. (2019), although with a small offset in the case of the latter due to the different high-z samples used. However, it shows a small discrepancy at M₂₀₀ ≳ 10^11.5 M_⊙ (a tiny part of the FL sample at z = 5.25) with the models of Moster et al. (2018). Some modern cosmological simulations with phenomenological feedback models (Pillepich et al. 2018; Davé et al. 2019) assume a decrease in the mass of galactic outflows with increasing galaxy mass to correct this discrepancy. However, these models are calibrated using observations at lower redshifts and may fail at earlier times, when gas densities were much higher. The FL suite, instead, predicts a higher density of massive galaxies (M^★ > 10⁹ M_⊙) resulting in an improved agreement with JWST data in the range z = 6 − 13, compared to the predictions of other models (Ceverino et al. 2024). On the other hand (bottom panel), the relation is in good agreement (except for an offset similar to that described above) with the results of the Cosmos-Web survey, obtained with both abundance-matching (Shuntov et al. 2025) and HOD (Paquereau et al. 2025) techniques. The tight relation observed in Fig. 2 implies that the FL virial masses, and therefore the virial radii, are typical of their stellar masses. We conclude that the compactness observed in a subpopulation of FL galaxies is due to a stellar distribution anomalously concentrated given its mass³. This sub-population must therefore arise from baryonic mechanisms.

Fig. 2. Relation between the stellar mass of galaxies and the virial mass of their haloes, compared to those predicted by semi-empirical models (top panel) and the Cosmos-Web survey subproducts (bottom panel). The shaded regions represent the dispersion of the models.

To shed light onto the origin of the compact subpopulation, we explored the evolution of the sizes of the 169 individual galaxies in our catalogue. Fig. 3 shows the four types of behaviour found: (a) galaxies that always expand, (b) objects showing an expansion at early times followed by a compaction at late times, (c) systems compacting first and expanding later, and (d) a combination of (b) and (c). We interpret this variety of behaviours as different stages of a single process in which galaxies first expand, then contract, and finally re-expand (hereafter ECE process). The compact subpopulation comprises, then, those galaxies going through the end of the contraction stage and the beginning of the following re-expansion.

Fig. 3. Different behaviours observed in the evolution of the sizes of individual galaxies: always expanding objects (top row), those reverting expansion into compaction (second row), systems reverting compaction into expansion (third row), and galaxies showing the complete ECE process (bottom row). Dashed and dash-dotted lines mark the masses at which galaxies attain their maximum and minimum sizes, respectively. Each panel represents a galaxy, labelled by its simulation ID.

The distribution of the masses $M_{\max }^{★}$ at which galaxies attain their peak sizes (Fig. 4) reveals that compaction predominantly occurs around M^★ ∼ 10^8 − 9 M_⊙. We define the turn-on mass, $M_{\mathrm{on}}^{★}$, as their geometric mean,

Fig. 4. Distribution of the masses at which galaxies attain their peak (left panel) and minimum (right panel) sizes. The dashed and dash-dotted lines represent the geometric mean values $M_{\mathrm{on}}^{★}$ and $M_{\mathrm{off}}^{★}$, respectively.

$$\begin{array}{c}\hfill log(M_{\mathrm{on}}^{★}/{\mathrm{M}}_{\odot })=⟨log(M_{\max }^{★}/{\mathrm{M}}_{\odot })⟩=8.53.\end{array}$$(1)

This represents the mass scale at which galaxies revert their initial expansion and begin to contract. The masses $M_{\min }^{★}$ at which galaxies begin re-expansion are concentrated around M^★ ∼ 10^9 − 10 M_⊙ (Fig. 4). We define the mean of this distribution as the turn-off mass.

$$\begin{array}{c}\hfill log(M_{\mathrm{off}}^{★}/{\mathrm{M}}_{\odot })=⟨log(M_{\min }^{★}/{\mathrm{M}}_{\odot })⟩=9.57.\end{array}$$(2)

We remark that the different phases of the ECE process occur at specific mass scales, indicating that galaxy stellar mass is its main driver.

To better analyse the development of this process, we divided our sample of all 7605 snapshots into four SSs:

• SS1: includes galaxies that always expand. These serve as a reference sample, representing systems not affected by the process within the redshift interval studied.

• SS2: comprises galaxies during the initial expansion phase, prior to compaction.

• SS3: includes galaxies undergoing the compaction phase, characterised by size reduction.

• SS4: consists of galaxies in the re-expansion phase, following the compaction event.

Note that, by construction, all galaxies except those in SS1 contribute to more than one SS; we include in each SS only the redshifts at which the galaxy displays the corresponding behaviour.

Fig. 5 shows the stellar size–mass relation of our sample, for redshifts z ∈ [5.25, 9]. Two different behaviours are already apparent at all redshifts: a tight ascending branch in which galaxy sizes increase with their masses (hereafter normal branch, because it follows the same trend observed at low redshifts, Shen et al. 2003; van der Wel et al. 2014; Mowla et al. 2019b), and a descending one with the opposite behaviour and a dispersion that increases with decreasing redshift. The break separating the two branches occurs at $M^{★}\approx M_{\mathrm{on}}^{★}$, suggesting that the change in the trend is connected to the beginning of the compaction phase of galaxies. Indeed, galaxies in SS1 and SS2 generally lie at masses lower than $M_{\mathrm{on}}^{★}$, whereas the opposite is true for SS3 and SS4, clearly showing once again that this mass scale is a global proxy for the onset of compaction. Above $M_{\mathrm{on}}^{★}$, only 18% of galaxies are expanding, whereas 58% are contracting and the remaining 24% have already contracted and are experiencing re-expansion. A second ascending branch due to the re-expansion of galaxies is not observed because the scatter in $M_{\min }^{★}$ mixes compacting and re-expanding galaxies at a given mass, producing a dispersion in the descending branch instead. A wider redshift range, allowing for the exploration of the whole re-expansion phase, would be needed to determine if a second ascending branch does exist.

Fig. 5. Dependence of the stellar size on stellar mass for different redshifts and galaxy SSs (indicated with different colours, SS1: red, SS2: yellow, SS3: green, SS4: blue). The dashed vertical line indicates the value of the turn-on stellar mass, $M_{\mathrm{on}}^{★}={10}^{8.53}{\mathrm{M}}_{\odot }$.

These trends are clear also in the relation between the stellar-to-halo size ratio, ξ^★, and the stellar mass (Fig. 6). This relation has a negative slope at all masses, originated in the already discussed faster growth of haloes with respect to the stellar distribution. A break separates the low- and high-mass branches at $M_{\mathrm{on}}^{★}$, the latter being steeper than the former. All galaxies with masses below $M_{\mathrm{on}}^{★}$ belong to the regular subpopulation (log ξ^★ > −1.6); compact galaxies are found only in the high-mass branch.

Fig. 6. Dependence of stellar-to-halo size ratio on stellar mass for different redshifts and galaxy SSs. References are as in Fig. 5. The dotted horizontal lines separate regular from compact galaxies.

In order to explore the physical nature of the compaction mechanism and what conditions trigger it, in Fig. 7 we plot the dependence of different baryon properties on stellar mass. The fraction of baryons and cold gas (i.e. gas at temperature T < 10⁴ K) within radius r are defined as f_b(r) = M_b(r)/M_t(r) and f_cold(r) = M_cold(r)/M_g(r), respectively. Here M_b(r), M_t(r), M_g(r), and M_cold(r) are, respectively, the total, baryonic, gas, and cold gas masses enclosed inside r. The mean temperature of cold gas, T_cold(r), was computed as the density-weighted average within the same volume, whereas its mean density is n_cold(r) = 3M_cold(r)/4πr³m_H, with m_H the mass of the hydrogen atom.

Fig. 7. Dependence of different properties of baryons on the stellar mass. Top panels: Fraction of baryons (left) and cold gas (right) within $r_{\mathrm{hm}}^{★}$. Bottom panels: Density-weighted mean gas temperature (left) and mean density (right) within $r_{\mathrm{hm}}^{★}$. The grey shades represent the density of galaxies in the plane defined by stellar mass and the corresponding property, computed from the data using a Gaussian smoothing kernel (normalization is different for each panel). Colour contours enclose 80% of the systems in different SSs (SS1: red, SS2: yellow, SS3: green, SS4: blue). The dashed and dash-dotted vertical lines mark $M_{\mathrm{on}}^{★}$ and $M_{\mathrm{off}}^{★}$, respectively, whereas the dotted horizontal line represent the density threshold for star formation.

When evaluating these properties within $r_{\mathrm{hm}}^{★}$, there are clear differences between galaxies of different SSs. The sequence SS1–SS2–SS3–SS4 shows progressively higher baryon fractions from $f_{\mathrm{b}}(r_{\mathrm{hm}}^{★})\sim 0.25$ up to almost unity (top left panel of Fig. 7), with a large jump at the turn-on mass scale. This suggests that the ECE process is related to the transition of the central regions of galaxies from dark-matter-dominated to baryon-dominated systems. It might be argued that the change in the baryon fraction could be interpreted as a consequence of the inward displacement of $r_{\mathrm{hm}}^{★}$, because central regions contain more baryons. However, an increase in $f_{\mathrm{b}}(r_{\mathrm{hm}}^{★})$ with M^★ is already seen in the reference SS1, which comprises galaxies that never contract.

Our results point therefore to the existence of an infalling flow of gas reaching the central regions, which modifies the ratio of baryons to dark matter, in agreement with previous works (Zolotov et al. 2015). Under this interpretation (for which we provide further evidence below), Fig. 7 (top right panel) shows also that the fraction of the cold gas phase increases from $f_{\mathrm{cold}}(r_{\mathrm{hm}}^{★})\sim 0.9$ up to ∼1. This happens because the infalling flow comprises mainly cold gas and/or the enhanced density accelerates cooling. Moreover, during compaction (SS3) and the following re-expansion (SS4), almost all gas in the central regions of galaxies is cold. This gas has reached conditions for star formation (temperature T < 10⁴ K, density n_cold > 1 cm⁻³, bottom panels of Fig. 7) in almost all systems except for the less massive ones in SS1. Indeed, this medium becomes cooler and denser as the ECE process develops, which reduces the local free-fall time, and thus shortens the star formation timescale. As a result, SFR density is strongly enhanced.

Our results point therefore to a wet compaction scenario (Dekel & Burkert 2014; Zolotov et al. 2015), characterised by the infall of large amounts of cold gas from the outer regions of galaxies, which compresses and cools up to the point that a starburst develops in the innermost regions. If the starburst is strong enough and occurs near the galaxy centre, it may reverse the normal inside-out star formation, turning the expansion of the stellar distribution into compaction. However, in previous works compaction typically occurs at lower redshifts (z ∼ 2 − 4) and a higher critical stellar mass (∼10¹⁰ M_⊙; Dekel & Burkert 2014; Zolotov et al. 2015; Tacchella et al. 2016). While driven by the same underlying mechanisms (self-reinforced cold gas inflows and central starburst activity) the ECE process identified in our simulations takes place earlier in the history of the Universe (9 > z > 5) and at a lower mass scale. It may indeed be responsible for the formation of spheroidal structures that would become the bulges of low-redshift galaxies (cf. Roper et al. 2022). After compaction, conditions may again favour star formation over a more extended region, resuming the inside-out growth of galaxies. Given that our ECE process requires that the stellar mass in the central regions exceeds the dark mass, and that the former is directly related to the total stellar mass, it is hardly conceivable that the differences with previous works, especially the critical masses, are due to different physics or calibration issues. Most probably we are seeing different stages in the evolution of the galaxies at which the wet compaction mechanism plays an important role.

Three issues must be explored to underpin the scenario described above. First, the time and spatial evolution of the process must agree with the described scenario. Second, the reason behind the efficient gas cooling that makes the starburst massive must be found. Third, the mechanism driving the strong gas infall must be pinpointed. Fig. 8 shows the fraction of baryons, $f_{\mathrm{b}}^{\mathrm{shell}}$, and cold gas, $f_{\mathrm{cold}}^{\mathrm{shell}}$, in a shell with inner and outer radii $r_{\mathrm{hm}}^{★}$ and 0.15 r₂₀₀, respectively, compared to the same fractions within $r_{\mathrm{hm}}^{★}$. The outer shell radius is well beyond $r_{\mathrm{hm}}^{★}\lesssim 0.1r_{200}$ for all galaxies (Fig. 1) and roughly marks the outer boundary of the stellar distribution. The ratios are near unity in SS1, demonstrating the existence of similar conditions in the central and outer regions of galaxies. They progressively drop by a significant amount in SS2–SS3–SS4, indicating that the key changes in baryonic properties responsible for driving the ECE process take place exclusively in the central regions.

Fig. 8. Ratio of the fraction of baryons (left panel) and cold gas (right panel) in a shell between $r_{\mathrm{hm}}^{★}$ and 0.15 r200 to those within $r_{\mathrm{hm}}^{★}$, as a function of stellar mass. References are as in Fig. 7. Dotted horizontal lines mark unit ratios.

The picture described above is supported also by the analysis of the evolution of the radial distribution of baryons in individual galaxies. Fig. 9 (top left panel) shows the radial flux of cold gas, Ṁ_infall, (defined as positive when gas flows inwards) as a function of radius and stellar mass, for a typical galaxy during the SS2–SS3–SS4 phases. It confirms the development of a strong infalling flow when the system reaches a stellar mass of the order of $M_{\mathrm{on}}^{★}$. The flux is significant only well within $r_{\mathrm{hm}}^{★}$ during the whole evolution, thus discarding that the changes we are seeing in baryon properties are merely the consequence of looking closer to the centre as $r_{\mathrm{hm}}^{★}$ shrinks. This flow carries cold gas to the inner kiloparsec, increasing the density of star-forming gas (ρ_g, SF; Fig. 9, top right panel), and thus triggering a localised star formation episode that leads to a central cusp in the SFR density profile (ρ_SFR; Fig. 9, bottom left panel). The gas densities ensure a high SFR via the Schmidt-Kennicutt law, making newborn stars dominate the stellar population over old ones. This produces a more compact stellar density distribution (ρ^★; Fig. 9, bottom right panel) and therefore a lower $r_{\mathrm{hm}}^{★}$.

Fig. 9. Radial profiles of the infalling mass flux (top left panel), and the densities of star-forming gas (top right panel), SFR (bottom left panel), and stars (bottom right panel), for an individual galaxy across its evolution. Colours refer to the evolving stellar mass of the galaxy. The inset shows the evolution of the galaxy size.

The aforementioned scenario explains the compaction phase of the process as originated in a self-reinforced infall episode that drives a sudden change of the prevailing gas conditions in the inner regions of galaxies. The recovery of an expansion phase, instead, is not led by a subsequent reverse change in these conditions. Galaxies in SS4 follow the same trends of those in SS2–SS3, and indeed represent extreme values with respect to the latter (Fig. 7). The profiles in Fig. 9 show not only that the infall persists at late times, but also that the flux reaches higher values. The cold gas, SFR, and stellar densities continue to increase at the centre, but develop growing tails up to larger radii. This suggests that the expansion is due to the enlargement of the region in which cold gas reaches star formation conditions. This re-expansion is different from that proposed by previous works (Dekel & Burkert 2014; Zolotov et al. 2015; Tacchella et al. 2016). Although these authors find a similar mass scale for re-expansion, they argue that this phase is driven by central gas depletion, which leads to the quenching of star formation in the core. Stars are then formed in a ring fuelled by fresh gas accreted at larger radii, resuming the inside-out growth of the galaxy. The fact that our galaxies reach the critical mass scale at earlier times suggests that they are richer in gas, which makes depletion less probable.

Fig. 10 shows the evolution of the gas depletion and replenishment times in the central regions, $t_{\mathrm{dep}}=M_{\mathrm{g}}(r_{\mathrm{hm}}^{★})/{\dot{M}}_{\mathrm{SFR}}(r_{\mathrm{hm}}^{★})$ and $t_{\mathrm{rep}}=M_{\mathrm{g}}(r_{\mathrm{hm}}^{★})/{\dot{M}}_{\mathrm{infall}}(r_{\mathrm{hm}}^{★})$, respectively, where Ṁ_SFR(r) is the total SFR rate within radius r. Expanding galaxies (SS1–SS2) display progressively lower depletion times due to the decrease of their gas reservoir, in line with the findings of Ceverino et al. (2018). However, t_dep falls barely below the timespan of the simulations at the end of this stage, indicating that this reservoir is never completely exhausted. The decrease in t_dep stops as soon as compaction starts, driven by the fresh gas infall. During the compaction and re-expansion stages (SS3–SS4) the depletion time remains almost constant, whereas the replenishment time decreases as infall becomes increasingly relevant. The two times become almost equal at the latest stage (SS4), indicating that the gas in the central region is replaced by fresh material from outside at a similar rate that it is used in star formation, and therefore never completely depleted. The existence of outflows may enhance depletion but would hardly modify this result, as outflow rates are usually similar to the SFR (Zolotov et al. 2015).

Fig. 10. Depletion time (left panel) and the ratio of replenishment-to-depletion times (right panel) within $r_{\mathrm{hm}}^{★}$, as a function of stellar mass. References are as in Fig. 7. The dotted horizontal line in the left panel shows the timespan of the simulations, whereas that in the right panel marks the unit ratio.

Fig. 11 provides clear evidence regarding the conditions that make possible the strong starburst found in our previous results. For the same typical galaxy of Fig. 9, the increase in the mass of the central regions driven by the infall of baryons (Fig. 11, top left panel) causes a decrease in the free-fall time, t_ff, in the same region (Fig. 11, top right panel). This enhances the infall itself, creating a self-reinforced flow of gas that increases the density in the central regions by several orders of magnitude. This gas cools in a timescale of

Fig. 11. Radial profiles of the fraction of baryons (top left panel), the free-fall (top right panel) and cooling times (bottom left panel), and the ratio tcool/tff (bottom right panel), for the same galaxy of Fig. 9. The dashed horizontal line marks tcool = tff.

$$\begin{array}{c}\hfill t_{\mathrm{cool}}=\frac{3.45k_{\mathrm{B}}T}{n_{\mathrm{g}}\mathrm{\Lambda }(n_{\mathrm{g}},T,Z)},\end{array}$$(3)

where k_B is Boltzmann constant, n_g and Z are the number density and metallicity of the gas, respectively, and Λ is the cooling rate per unit density. The cooling rate shows a strong dependence of n_g (Wiersma et al. 2009), ensuring that in the central regions the cooling time is short (Fig. 11, bottom panels) and almost all gas remains cold. The increase in density produced by the infall thus causes a massive starburst that makes the stellar distribution more compact.

Fig. 12 provides further evidence for our scenario by displaying the dependence on stellar mass of the density-weighted average radii of the cold gas and SFR distributions (r_cold and r_SFR, right and left panels, respectively). These radii estimate the size of the regions where the cold gas lies and the stars are being formed, respectively. They were computed as

Fig. 12. Dependence on the stellar mass of the density-weighted average radii of the distribution of cold gas (top left panel) and SFR (top right panel). References are as in Fig. 7.

$$\begin{array}{c}\hfill r_{\mathrm{cold}}={⟨r⟩}_{{\rho }_{\mathrm{cold}}}=\frac{{\int }_0^{0.15r_{200}}r{\rho }_{\mathrm{cold}}(r)\mathrm{d}r}{{\int }_0^{0.15r_{200}}{\rho }_{\mathrm{cold}}(r)\mathrm{d}r},\end{array}$$(4)

and

$$\begin{array}{c}\hfill r_{\mathrm{SFR}}={⟨r⟩}_{{\rho }_{\mathrm{SFR}}}=\frac{{\int }_0^{0.15r_{200}}r{\rho }_{\mathrm{SFR}}(r)\mathrm{d}r}{{\int }_0^{0.15r_{200}}{\rho }_{\mathrm{SFR}}(r)\mathrm{d}r}·\end{array}$$(5)

The upper limit of integration was set to 0.15 r₂₀₀, which approximately corresponds to the outer boundary of the stellar distribution. Both r_cold and r_SFR are close to $r_{\mathrm{hm}}^{★}$ for galaxies in SS1, indicating that the star formation is widely distributed, which is consistent with the slow expansion of these galaxies. These radii fall well below $r_{\mathrm{hm}}^{★}$ in SS2 and SS3, indicating that the star formation concentrates well inside $r_{\mathrm{hm}}^{★}$, as proposed. Moreover, galaxies in SS2 are still expanding; the decrease in the half-mass radius lags behind those of the other radii, which is consistent with the latter decrease being the cause of the former. Fig. 12 also shows that the cold gas and SFR distributions of SS4 galaxies are more extended with respect to the half-mass radius than those of SS3, indicating that star formation conditions proceed in SS4 into a broader region, which causes the recovery of an expansion phase.

Finally, the mechanism behind the self-reinforced infall–cooling–star formation episode can be unveiled from Fig. 13. The infalling flux, Ṁ_infall, at $r_{\mathrm{hm}}^{★}$ is roughly dependent on the square of the total mass, M_tot, within this radius, which is reminiscent of Bondi–Hoyle accretion (Bondi & Hoyle 1944). Indeed, a simple model (dashed lines in Fig. 13) constructed by assuming that the flow proceeds in this regime provides a reasonably good agreement with the data at all redshifts. In this model,

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{infall}}=\frac{\pi \rho {(1+z)}^3{G}^2{M}_{\mathrm{tot}}^2}{{(k_{\mathrm{B}}T{m}_{\mathrm{H}}^{-1})}^{3/2}},\end{array}$$(6)

where G is the gravitational constant and m_H is the mass of the hydrogen atom. We assume that gas reaches $r_{\mathrm{hm}}^{★}$ from outside at a typical temperature T = 10 000 K and density ρ = f_b, typΩ_bρ_crit, where ρ_crit, 0 = 9.47 × 10⁻²⁷ kg m⁻³ is the critical density of the Universe at redshift z = 0, and f_b, typ = 0.05 was taken as the mean baryon fraction of SS1 galaxies at 0.15 r₂₀₀. While the galaxies remain dark-matter dominated, the infalling flux depends mainly on the dark matter mass in the inner regions. As soon as baryons dominate, the flux couples to the baryonic mass provided by itself, and the dependence of Bondi–Hoyle accretion on M_tot² generates the self-reinforced process observed. This change of behaviour happens roughly at the turn-on mass scale.

Fig. 13. Flux Ṁinfall of infalling gas at $r_{\mathrm{hm}}^{★}$ as a function of total mass, Mtot, within this radius, for different redshift intervals and galaxy SSs (indicated with different colours, SS1: red, SS2: yellow, SS3: green, SS4: blue). Dashed lines represent a semiempirical model for the data, describing the infalling flux of gas in the Bondi–Hoyle regime.

This accretion mechanism sheds light also on the development of the re-expansion phase. Fig. 14 shows the change of the mean baryon fraction, f_b(r), within spheres of progressively larger radii r. It is clear that the transition of central regions from dark-matter- to baryon-dominated proceeds in an inside-out way. The process starts near M_on only at the innermost regions, which increase their f_b and develop a sizeable inflow of gas that leads to a strong, concentrated starburst. This causes the stellar half-mass radius to shrink, as discussed above. As baryon dominance spreads outwards the self-reinforced accretion flow switches on at progressively larger radii, leading to a wider starburst and a less concentrated distribution of newly born stars. This stops and eventually reverses the decrease in $r_{\mathrm{hm}}^{★}$, driving the observed re-expansion.

Fig. 14. Fraction of baryons within spheres of different radii r (indicated by different colours) as a function of stellar mass. Solid lines display the median over all galaxies, whereas shaded regions represent the dispersion of each distribution.

In summary, the three issues raised above have been settled. The time and spatial evolution of our simulated galaxies strongly support the wet compaction scenario, with a final re-expansion driven by the enlargement of the star-forming region instead of central gas depletion. Gas infall into the central regions of galaxies is the main driver of the process; the dependence of the infalling flux on the square of the mass of the accreting region is its key. At low stellar masses and baryon fractions, the flux depends mainly on the dominant dark-matter mass, and therefore proceeds at a slow rate. As soon as this slow infall begins to transform the central region of a galaxy into a baryon-dominated system, the dominant role is taken by the baryon mass. Infall becomes then faster through self-reinforcement, as it increases the same baryon mass that drives it. The result is a strong enhancement of the gas density in the central region, a large drop in the gas cooling time, and the consequent increase in the SFR. These produce a strong, concentrated starburst that makes the stellar distribution more compact.

4. Discussion and conclusions

Recent JWST observations have challenged our understanding of galaxy formation during the Cosmic Dawn. Specifically, the physical mechanisms creating compact galaxies with sizes of the order of 100 pc, and those leading to an inverted stellar mass–size relation at this epoch, are still a matter of debate. Aiming at unveiling those mechanisms, in this paper we have analysed the evolution of the stellar distribution of 169 galaxies from the FL zoom-in simulation suite, which enables a statistically robust analysis of these early evolutionary processes.

We found that, underneath the bulk evolution due to hierarchical growth, a population of compact galaxies with sizes of the order of ∼100 pc develops. This population is a result of a process that stops the normal inside-out growth of galaxies, suddenly shrinking their stellar distribution to sizes comparable to those observed by JWST (Baggen et al. 2023; Baker et al. 2025; Tacchella 2025; Morishita et al. 2024). Some of the observations are indeed compatible with our results in both size and mass (see Fig. 4 of Morishita et al. 2024). After this rapid compaction, the stellar mass distribution resumes its inside-out development. Remarkably, both the onset of compaction and the resumption of expansion occur near specific stellar mass scales, $M_{\mathrm{on}}^{★}={10}^{8.53}{\mathrm{M}}_{\odot }$ and $M_{\mathrm{off}}^{★}={10}^{9.57}{\mathrm{M}}_{\odot }$, respectively, suggesting that stellar mass is the main driver of the process. This happens because compaction begins when the stellar mass in the central regions, directly related to the total stellar mass, matches the dark mass in the same regions. Variations in the critical masses may therefore be due to different individual star formation histories and spatial distributions produced by, for example, the environment or the merger history. Indeed, the mechanism operates in all galaxies that reach the corresponding mass scales during the analysed redshift interval. Galaxies in which we do not observe it have masses below $M_{\mathrm{on}}^{★}$ at the lowest redshift, whereas galaxies with masses above $M_{\mathrm{off}}^{★}$ at the highest redshift are all compact. This leads us to propose that the process is universal; which may be investigated further by enlarging our redshift interval and simulation box. We are working on this issue at present.

The described process naturally leads to an inversion of the size–mass relation in the stellar mass range 10^8.5 − 10^9.5 M_⊙. This result is consistent with the observations of Baggen et al. (2023), and with recent numerical works based on simulations such as BLUETIDES at z = 7 (Roper et al. 2022; Marshall et al. 2022). However, it is at odds with the conclusion of Morishita et al. (2024) who report a positive slope in the whole mass range, and also with other zoom-in simulations such as THESAN-ZOOM (McClymont et al. 2025) and FIREBOX (Feldmann et al. 2023; Rohr et al. 2022), whose results are consistent with such trend. The inversion in the mass-size relation has been previously reported in numerical simulations by Zolotov et al. (2015) and Lapiner et al. (2023). Both studies observed that, for certain mass ranges, the size exhibits a decrease with mass followed by the typical increase. In the case of Lapiner et al. (2023), the mass ranges corresponding to this behaviour are similar to those found in our study; in other cases it is slightly higher (10¹⁰ M_⊙). A notable difference, however, lies in the redshift coverage: Lapiner et al. (2023) explores the range from z ≈ 1 to 5, while Zolotov et al. (2015) focuses on z ≈ 2 to 4. In this regard, our work provides new insights by extending the analysis to higher redshifts, from z ≈ 5 to 9, with good resolution. This extended range allows us to capture the onset of the compaction process, and therefore the full ECE behaviour – our sample is capable of tracing the entire cycle of size evolution. The agreement of the results throughout the whole interval 1 ≲ z ≲ 9 supports our claim that this process is universal.

Regarding the physics, our results point to a wet compaction scenario similar to those described in previous works (Dekel & Burkert 2014; Zolotov et al. 2015). The infall of gas from the outer skirts of a galaxy to its innermost (< 1 kpc) region increases the density of the latter and changes its state from dark matter- to baryon-dominated. It also decreases the cooling time of the gas, enhancing the fraction of cold gas, and boosts star formation. The process becomes self-reinforced because the infalling flux increases with the baryon mass in the inner regions. The result is a strong central starburst, with newborn stars dominating the whole stellar distribution and making it appear more compact.

Although the development of gas cooling and star formation conditions proceeds inside-out, the starburst is strong and concentrated enough to produce a rapid and steady decrease in the stellar half-mass radius. The subsequent re-expansion phase is driven by the broadening of the region in which baryons dominate, a natural consequence of gas flowing inwards, and the resulting enhancement of accretion at progressively larger radii. Therefore, cold gas can reach star formation conditions farther out from the centre, leading to a wider stellar distribution and a larger half-mass radius. This is a major difference with previous works, in which it is due to the depletion of gas and the corresponding quenching of star formation in the central regions, which leads to a star-forming ring that grows outwards. We argue that the galaxies at low redshift, such as those studied in previous works, have low gas densities that facilitate quenching. Our high-redshift galaxies are gas-rich, which makes it difficult to exhaust their gas; it continues to form stars while more gas reaches the inner regions, expanding the star-forming region. Given that re-expansion is attained at the same mass scale but different redshifts, the development or quenched or normal re-expansion might depend on the environment. This issue deserves further exploration.

More recently, using the THESAN zoom-in simulations, McClymont et al. (2025) find that the evolution of galaxy sizes is closely linked to bursty star formation at high redshift. In their simulations, galaxies located above the star-forming main sequence tend to exhibit centrally concentrated star formation, which drives morphological compaction episodes. These central starbursts subsequently quench in an inside-out manner, resulting in cyclical alternations between compaction and expansion that cause galaxies to oscillate around an overall normal size–mass relation, instead of a producing a steady compaction followed by an expansion. McClymont et al. (2025) claim that the reason behind the discrepancy between their results and those of other simulations is their greater mass resolution that produces a set of small starbursts instead of a single, massive one. It may be argued, however, that the overall increase in the stellar mass of galaxies is the same in all simulations; therefore the different behaviour of $r_{\mathrm{hm}}^{★}$ is not due to mass resolution but to the varying spatial distribution of the starbursts with time. Fig. 4 of McClymont et al. (2025) shows starbursts changing the half-mass radius by ∼ − 0.5 dex, followed by a rapid recovery of its previous value. This implies a large mass of newborn stars is initially created at the centre, but also that a comparable mass is produced outside $r_{\mathrm{hm}}^{★}$ afterwards. Therefore, quenching must proceed rapidly outwards and involve a high mass and a large spatial scale, both well above resolution. Moreover, their half-mass radii are $r_{\mathrm{hm}}^{★}\gtrsim 1$ kpc, indicating a more extended star-formation region than in both our simulation and observed galaxies. The discrepancy may therefore probably be related to different physical modelling. Indeed, our resolution is comparable to the lowest one in the THESAN simulations. We note that extreme resolution values of a few tens of parsecs and hundreds of solar masses require a thorough evaluation of sub-grid processes, especially star formation, as some base assumptions such as the Kennicutt–Schmidt law are not tested down to such scales.

In principle, the gas flows that lead to compaction can be affected by stellar feedback. Stronger stellar feedback and larger velocity dispersion in the galaxies may disrupt the compaction described here until larger masses can be reached. However, different simulations with very different feedback models all describe a qualitatively similar central or bulge growth (Ceverino et al. 2010; Genel et al. 2012; Mandelker et al. 2017; Oklopčić et al. 2017). It seems that dense gas inflows at galaxy-scale are very resilient to the injection of energy and/or momentum. In particular, Ceverino et al. (2023) discuss the effect of two different feedback models on the compaction at relatively low redshifts (z = 1 − 3). The evolution of compaction and early quenching is similar in both cases, although stellar masses (and central densities) are a factor of two lower in the strong feedback case. On the other hand, the addition of AGN feedback may affect the growth of massive galaxies, especially during their last re-expansion phase (Lapiner et al. 2023).

On the other hand, a careful comparison between observations and simulations involves several critical details that need to be addressed. In this context, it is essential to fully consider the observational procedures. We are currently addressing this issue in a forthcoming paper. It is not clear if an inverted size–mass relation should be directly observable due to the spread in the masses at which galaxies start and finish compaction. Drivers responsible for this spread should be explored in order to identify galaxy SSs that may cleanly show the inversion, before a robust conclusion can be drawn from the observed size–mass relation. Indeed, although not a physical driver, our results show that redshift may indirectly play this role: the higher the redshift, the cleaner the inverted size–mass relation appears in our simulations. A large sample of galaxies in the aforementioned mass range and beyond z ∼ 8 may provide clues on this issue.

Another issue is the transformation of the physical size–mass relation into an observable one. Observers tackle the inverse problem through broadband SED fitting (e.g. Roper et al. 2022; Morishita et al. 2024), parameterising their models with the stellar mass and star formation history, and including the emission of ionised gas and dust absorption. The masses obtained are therefore subject to uncertainties and degeneration present in stellar evolutionary models, radiation transfer physics, and so on. On the other hand, the sizes are measured from projected luminosity profiles, which are not easily deconvolved into 3D light distributions, and are subject to observational biases. A different (direct) approach is to construct galaxy images from the simulations, and compare directly the size–mass relations obtained from them. Biases can be included in the construction, and their effects on the relation explored. We will be presenting the results of such an analysis in a forthcoming paper.

Finally, an important characteristic of the proposed self-reinforced gas infall to the centre of galaxies is that it is self-triggered. We have shown that the accretion of this gas proceeds in the Bondi & Hoyle (1944) regime, the flow being dependent of the square of the accreting mass. As far as the accreted mass overcomes that originally present in this region, the infall is dominated by the former and the process becomes self-reinforced. Indeed, we observe this change to occur because the inflowing matter is baryon-dominated, whereas the original mass in the central regions is dominated by dark matter. The switching between both regimes (from the prevalence of dark matter to that of baryons) is clearly seen in our simulations and supports our analysis of the physical mechanism behind compaction. It also explains why mass is its main driver. The Bondi & Hoyle (1944) regime describes spherical accretion, in which angular momentum of the accreted material can be neglected. The ECE process may therefore constitute a continuation of the outer angular momentum loss process highlighted by Lapiner et al. (2023). It may also be responsible for the formation of spheroidal structures with no or little net angular momentum, which would become the bulges of low-redshift galaxies, as has been proposed by Roper et al. (2022). As the re-expansion phase begins and star formation moves to larger radii, angular momentum should play a more relevant role, changing or stopping the infall regime, and favouring the emergence of rotation-dominated structures such as discs. An investigation of the morphology and dynamics of the structures created during the ECE process is under way, and will be presented elsewhere.

In summary, our main conclusions are the following:

• We have described the complete compaction process of galaxies since its onset to the late re-expansion, providing evidence for the wet compaction scenario put forth by previous works.

• Together with previous results, our findings suggest that this scenario is universal, operating at all redshifts.

• The re-expansion ending compaction proceeds without central quenching at high redshift, contrary to its low-redshift analogue. This may be due to the higher gas density of galaxies at earlier epochs.

• We provided a physical explanation for the onset of wet compaction and its mass scale, based on Bondi–Hoyle accretion. In addition, we introduced a new analytical expression, potentially valuable for use in semi-analytical models.

Data availability

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

Acknowledgments

The authors acknowledge insightful comments by the anonymous referee, which helped to greatly improve the manuscript. PC, SP, LJP, and LAB acknowledge partial support from ANPCyT through grant PICT 2020-00582. This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement No 734374 (LACEGAL). DC is supported by research grant PID2021-122603NB-C21 funded by the Ministerio de Ciencia, Innovación y Universidades (MICIU/FEDER) and the research grant CNS2024-154550 funded by MICIU/AEI/10.13039/501100011033. The authors gratefully acknowledge the Gauss Center for Supercomputing for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Center (Project ID: pr92za). The authors thankfully acknowledge the computer resources at MareNostrum and the technical support provided by the Barcelona Supercomputing Center (RES-AECT-2020-3-0019). We acknowledge the publicly available programming language python, including the numpy, astropy, matplotlib, scipy, and h5py packages.

References

All Figures

	Fig. 1. Distribution of the half-mass stellar radius (left panel) and the stellar-to-halo size ratio (right panel), for 5.25 < z < 9 (colour coded). The dashed vertical line indicates the median value of ξ★ over the full sample and redshift range (⟨ξ★⟩ = 0.07). A growing subpopulation of compact galaxies is seen; the dotted vertical line at log ξ★ = −1.6 roughly illustrates the separation.
In the text

	Fig. 2. Relation between the stellar mass of galaxies and the virial mass of their haloes, compared to those predicted by semi-empirical models (top panel) and the Cosmos-Web survey subproducts (bottom panel). The shaded regions represent the dispersion of the models.
In the text

Fig. 3. Different behaviours observed in the evolution of the sizes of individual galaxies: always expanding objects (top row), those reverting expansion into compaction (second row), systems reverting compaction into expansion (third row), and galaxies showing the complete ECE process (bottom row). Dashed and dash-dotted lines mark the masses at which galaxies attain their maximum and minimum sizes, respectively. Each panel represents a galaxy, labelled by its simulation ID.

In the text

	Fig. 4. Distribution of the masses at which galaxies attain their peak (left panel) and minimum (right panel) sizes. The dashed and dash-dotted lines represent the geometric mean values $M_{\mathrm{on}}^{★}$ and $M_{\mathrm{off}}^{★}$, respectively.
In the text

	Fig. 5. Dependence of the stellar size on stellar mass for different redshifts and galaxy SSs (indicated with different colours, SS1: red, SS2: yellow, SS3: green, SS4: blue). The dashed vertical line indicates the value of the turn-on stellar mass, $M_{\mathrm{on}}^{★}={10}^{8.53}{\mathrm{M}}_{\odot }$.
In the text

	Fig. 6. Dependence of stellar-to-halo size ratio on stellar mass for different redshifts and galaxy SSs. References are as in Fig. 5. The dotted horizontal lines separate regular from compact galaxies.
In the text

Fig. 7. Dependence of different properties of baryons on the stellar mass. Top panels: Fraction of baryons (left) and cold gas (right) within $r_{\mathrm{hm}}^{★}$. Bottom panels: Density-weighted mean gas temperature (left) and mean density (right) within $r_{\mathrm{hm}}^{★}$. The grey shades represent the density of galaxies in the plane defined by stellar mass and the corresponding property, computed from the data using a Gaussian smoothing kernel (normalization is different for each panel). Colour contours enclose 80% of the systems in different SSs (SS1: red, SS2: yellow, SS3: green, SS4: blue). The dashed and dash-dotted vertical lines mark $M_{\mathrm{on}}^{★}$ and $M_{\mathrm{off}}^{★}$, respectively, whereas the dotted horizontal line represent the density threshold for star formation.

In the text

	Fig. 8. Ratio of the fraction of baryons (left panel) and cold gas (right panel) in a shell between $r_{\mathrm{hm}}^{★}$ and 0.15 r200 to those within $r_{\mathrm{hm}}^{★}$, as a function of stellar mass. References are as in Fig. 7. Dotted horizontal lines mark unit ratios.
In the text

	Fig. 9. Radial profiles of the infalling mass flux (top left panel), and the densities of star-forming gas (top right panel), SFR (bottom left panel), and stars (bottom right panel), for an individual galaxy across its evolution. Colours refer to the evolving stellar mass of the galaxy. The inset shows the evolution of the galaxy size.
In the text

	Fig. 10. Depletion time (left panel) and the ratio of replenishment-to-depletion times (right panel) within $r_{\mathrm{hm}}^{★}$, as a function of stellar mass. References are as in Fig. 7. The dotted horizontal line in the left panel shows the timespan of the simulations, whereas that in the right panel marks the unit ratio.
In the text

	Fig. 11. Radial profiles of the fraction of baryons (top left panel), the free-fall (top right panel) and cooling times (bottom left panel), and the ratio tcool/tff (bottom right panel), for the same galaxy of Fig. 9. The dashed horizontal line marks tcool = tff.
In the text

	Fig. 12. Dependence on the stellar mass of the density-weighted average radii of the distribution of cold gas (top left panel) and SFR (top right panel). References are as in Fig. 7.
In the text

	Fig. 13. Flux Ṁinfall of infalling gas at $r_{\mathrm{hm}}^{★}$ as a function of total mass, Mtot, within this radius, for different redshift intervals and galaxy SSs (indicated with different colours, SS1: red, SS2: yellow, SS3: green, SS4: blue). Dashed lines represent a semiempirical model for the data, describing the infalling flux of gas in the Bondi–Hoyle regime.
In the text

	Fig. 14. Fraction of baryons within spheres of different radii r (indicated by different colours) as a function of stellar mass. Solid lines display the median over all galaxies, whereas shaded regions represent the dispersion of each distribution.
In the text