We investigate galaxy sizes at redshift inlineformulatexmath idTM0001 notationLaTeXzgtrsim 6texmathinlineformula with the cosmological radiationmagnetohydrodynamic simulation suite THESANHR. These simulations simultaneously capture reionization of the largescale intergalactic medium and resolved galaxy properties. The intrinsic sizes inlineformulatexmath idTM0002 notationLaTeXrast 12texmathinlineformula of simulated galaxies increase moderately with stellar mass at inlineformulatexmath idTM0003 notationLaTeXMast  lesssim 108 rm Modottexmathinlineformula and decrease fast at larger masses resulting in a hump feature at inlineformulatexmath idTM0004 notationLaTeXMast sim 108 rm Modottexmathinlineformula that is insensitive to redshift. Lowmass galaxies are in the initial phase of size growth and are better described by a spherical shell model with feedbackdriven outflows competing with the cold inflowing gas streams. In contrast massive galaxies fit better with the disc formation model. They generally experience a phase of rapid compaction and gas depletion likely driven by internal disc instability rather than external processes. We identify four compact quenched galaxies in the inlineformulatexmath idTM0005 notationLaTeX95.5 rm cMpc3texmathinlineformula volume of THESAN1 at inlineformulatexmath idTM0006 notationLaTeXzsimeq 6texmathinlineformula and their quenching follows reaching a characteristic stellar surface density akin to the massive compact galaxies at cosmic noon. Compared to observations we find that the median ultraviolet effective radius inlineformulatexmath idTM0007 notationLaTeXRrm UVrm efftexmathinlineformula of simulated galaxies is at least three times larger than the observed ones at inlineformulatexmath idTM0008 notationLaTeXMast lesssim 109 rm Modottexmathinlineformula or inlineformulatexmath idTM0009 notationLaTeXMrm UVgtrsim 20texmathinlineformula at inlineformulatexmath idTM0010 notationLaTeX6 lesssim z lesssim 10texmathinlineformula. The population of compact galaxies inlineformulatexmath idTM0011 notationLaTeXRrm UVrm efflesssim 300 rm pctexmathinlineformula galaxies at inlineformulatexmath idTM0012 notationLaTeXMast sim 108 rm Modottexmathinlineformula is missing in our simulations. This inconsistency persists across many other cosmological simulations with different galaxy formation models and demonstrates the potential of using galaxy morphology to constrain physics of galaxy formation at high redshifts.
['Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
The THESAN project galaxy sizes during the epoch of reionization
{'The /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc project: galaxy sizes during the epoch of reionization': 'Xuejian Shen, 1 Mark Vogelsberger, 1 , 2 Josh Borrow, 1 , 3 Yongao Hu, 1 Evan Erickson, 1 Rahul Kannan, 4 Aaron Smith, 5 Enrico Garaldi, 6 , 7 Lars Hernquist, 8 Takahiro Morishita, 9 Sandro Tacchella, 10 , 11 Oliver Zier, 1 Guochao Sun, 12 Anna-Christina Eilers 1 and Hui Wang 1  \n- 1 Department of Physics & Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA\n- 2 The NSF AI Institute for Artificial Intelligence and Fundamental Interactions, Massachusetts Institute of Technology, Cambridge, MA 02139, USA\n- 3 Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA\n- 4 Department of Physics and Astronomy, York University, 4700 Keele Street, Toronto, ON M3J 1P3, Canada\n- 5 Department of Physics, The University of Texas at Dallas, SCI10 800 West Campbell Road, Richardson, Texas 75080, USA\n- 6 Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy\n- 7 Max-Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany\n- 8 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA\n- 9 IPAC, California Institute of Technology, MC 314-6, 1200 E. California Boulevard, Pasadena, CA 91125, USA\n- 10 Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK\n- 11 Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK\n- 12 CIERA and Department of Physics and Astronomy, Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'We investigate galaxy sizes at redshift 𝑧 ≳ 6 with the cosmological radiation-magneto-hydrodynamic simulation suite /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc(/h.pc/r.pc). These simulations simultaneously capture reionization of the large-scale intergalactic medium and resolved galaxy properties. The intrinsic sizes ( 𝑟 ∗ 1 / 2 ) of simulated galaxies increase moderately with stellar mass at 𝑀 ∗ ≲ 10 8 M ⊙ and decrease fast at larger masses, resulting in a hump feature at 𝑀 ∗ ∼ 10 8 M ⊙ that is insensitive to redshift. Low-mass galaxies are in the initial phase of size growth and are better described by a spherical shell model with feedback-driven outflows competing with the cold inflowing gas streams. In contrast, massive galaxies fit better with the disk formation model. They generally experience a phase of rapid compaction and gas depletion, likely driven by internal disk instability rather than external processes. We identify four compact quenched galaxies in the ( 95 . 5 cMpc ) 3 volume of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 at 𝑧 ≃ 6 and their quenching follows reaching a characteristic stellar surface density akin to the massive compact galaxies at cosmic noon. Compared to observations, we find that the median UV effective radius ( 𝑅 UV eff ) of simulated galaxies is at least three times larger than the observed ones at 𝑀 ∗ ≲ 10 9 M ⊙ or 𝑀 UV ≳ -20 at 6 ≲ 𝑧 ≲ 10. The population of compact galaxies ( 𝑅 UV eff ≲ 300 pc) galaxies at 𝑀 ∗ ∼ 10 8 M ⊙ is missing in our simulations. This inconsistency persists across many other cosmological simulations with different galaxy formation models and demonstrates the potential of using galaxy morphology to constrain physics of galaxy formation at high redshifts.  \nKey words: methods: numerical - galaxies: evolution - galaxies: formation - galaxies: high-redshift - galaxies: structure', '1 INTRODUCTION': 'The size of the stellar distribution of galaxies and its evolution with cosmic time provide valuable insights about the formation history of galaxies and the relationship with their dark matter (DM) haloes (e.g. Mo et al. 1998; Kravtsov 2013; Somerville et al. 2018). At low redshifts, the two main categories of galaxies (star-forming and quenched) show rather different dependencies between galaxy size and stellar mass/luminosity (e.g. Shen et al. 2003; Kauffmann et al. 2003; Franx et al. 2008; van der Wel et al. 2014), which encode their different evolutionary and assembly histories. Galaxy sizes at lower redshifts and their complicated dependencies on various galaxy properties are well reproduced in cosmological hydrodynamic simulations (e.g. Genel et al. 2018; Pillepich et al. 2019; Popping et al. 2022).  \nHowever, at high redshifts, the classification of galaxy morphology and the measurement of galaxy sizes are challenging due to the limited angular resolution of instruments. Prior to the James WebbSpace Telescope (JWST), galaxy size measurements have been mainly driven by the Hubble Space Telescope (HST) with the Advanced Camera for Surveys and the Wide Field Camera 3/IR channel on board (e.g. Bruce et al. 2012; Mosleh et al. 2012; van der Wel et al. 2014; Morishita et al. 2014; Allen et al. 2017; Mosleh et al. 2020). This has been pushed for Lyman Break galaxies (LBGs) out to 𝑧 ∼ 8 with HST legacy data (e.g. Oesch et al. 2010; Ono et al. 2013; Shibuya et al. 2015; Bouwens et al. 2022). These studies have found stable slopes and scatters of size-mass/luminosity relation at 0 ≤ 𝑧 ≲ 8 but significantly decreasing average sizes and increased star-formation rate surface density towards high redshifts. These are  \nin good agreement with disk formation theories of star-forming/latetype galaxies (e.g. Mo et al. 1998; Dutton et al. 2007; Kravtsov 2013). At 𝑧 ∼ 3, there are no significant differences between the size of starforming versus quiescent galaxies at the massive end (e.g. Faisst et al. 2017; Hill et al. 2017; Mowla et al. 2019) when the morphological transitions are thought to present.  \nThe early results from the Early Science Release and Cycle-1 observations have already demonstrated the remarkable capabilities of JWST . With its infrared sensitivity and angular resolution, JWST revealed morphologies of early galaxies at 𝑧 ≳ 6 down to scales of ∼ 100 pc throughout rest-frame ultraviolet (UV) to optical wavelengths (e.g. Yang et al. 2022; Robertson et al. 2023b; Tacchella et al. 2023a; Treu et al. 2023; Huertas-Company et al. 2023; Ormerod et al. 2023). Notably, JWST has revealed a population of extremely compact galaxies at 𝑧 ≳ 6 (e.g. Baker et al. 2023; Baggen et al. 2023; Morishita et al. 2024). These galaxies are more compact than the extrapolations of low-redshift measurements or theoretical model predictions. The formation channel of compact galaxies at very early stages remains unexplored. Many of the compact galaxies show signatures of heavily obscured active galatic nuclei (AGN; e.g. Greene et al. 2023; Harikane et al. 2023; Labbe et al. 2023; Kocevski et al. 2023; Matthee et al. 2023; Maiolino et al. 2023; Kokorev et al. 2024) with abundance significantly exceeding the extrapolation of UV-selected AGN in the Hubble era (e.g. Kulkarni et al. 2019; Niida et al. 2020; Shen et al. 2020b). They could bias the constraints on galaxy morphology (e.g. Harikane et al. 2023; Tacchella et al. 2023b) but also pose interesting questions regarding the co-evolution of supermassive black holes (SMBHs) with their host galaxies.  \nOn the other hand, theoretical predictions on galaxy sizes have been made by several numerical simulations, which have the unique power to predict spatially resolved galaxy properties. For example, Maet al. (2018a) studied galaxy morphology and sizes at 𝑧 > 5 using the FIRE-2 simulations (Ma et al. 2018b) and focused on the lowmass end. Marshall et al. (2022) studied galaxy sizes at 7 < 𝑧 < 12 using mock images generated from the Bluetides simulations and highlighted the impact of dust. Roper et al. (2022) studied galaxy sizes at 𝑧 ≥ 5 using the FLARES simulations (Lovell et al. 2021) and explored physical mechanisms of compact galaxy formation (Roper et al. 2023). In general, these studies found a negative correlation between intrinsic sizes of galaxies with stellar mass that is contradictory to the observed size-mass/luminosity relation (e.g. Shibuya et al. 2015; Kawamata et al. 2018). Concentrated dust attenuation could play a key role in reversing the trend in massive luminous galaxies and reconciling theoretical predictions with observations (e.g. Marshall et al. 2022; Roper et al. 2022; Popping et al. 2022). However, it is worth noting that most of these simulations are calibrated based on low-redshift observables and are subject to great uncertainties when extrapolating to high redshifts. For instance, the inhomogeneity of the UV radiation background during the epoch of reionization can affect the star-formation and metal enrichment of low-mass galaxies in non-trivial ways (Rosdahl et al. 2018; Katz et al. 2020; Borrow et al. 2022) compared to models assuming spatially uniform radiation background. The implementation and calibration of different stellar feedback mechanisms are far from converging between theoretical frameworks and become more uncertain at high redshifts (e.g. Gnedin 2014; Pawlik et al. 2017; Pallottini et al. 2022; Dekel et al. 2023; Ferrara 2023).  \nIn this study, we aim to provide a thorough analysis of galaxy sizes at 𝑧 ≳ 6 over a large range of galaxy stellar mass from ∼ 10 6 . 5 to 10 10 . 5 M ⊙ using both the flagship /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations and its high-resolution variants. As a suite of radiation-hydrodynamic simulations, /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc is able to self-consistently model the reionization  \nprocess and the galaxies responsible for it with unprecedented physical fidelity. We will first focus on the physical mechanisms driving galaxy size evolution and understanding the causal connection between size and other galaxy properties. We will then compare results with the most recent JWST constraints and explore implications for physics models for galaxy formation at high redshifts. The paper is organized as follows. In Section 2, we introduce the simulation suite. In Section 3, we introduce the analysis methods, including galaxy identification, classification, post-processing, and definitions of various galaxy properties. In Section 4, we study the galaxy intrinsic size-mass relation and the correlation of size with other galaxy or environmental properties. In Section 5, we explore the physical mechanisms that drive galaxy size evolution, especially the compaction of massive galaxies. In Section 6, we compare the simulation results with observations and discuss various theoretical and observational uncertainties. In Section 7, we present the conclusions of the paper.', '2 SIMULATIONS': 'The /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc project (Kannan et al. 2022a; Garaldi et al. 2022; Smith et al. 2022) is a suite of radiation-magneto-hydrodynamic simulations utilizing the moving-mesh hydrodynamics code AREPO (Springel 2010; Weinberger et al. 2020). Gravity is solved using the hybrid Tree-PM method (Barnes & Hut 1986). The hydrodynamics is solved using the quasi-Lagrangian Godunov method (Godunov & Bohachevsky 1959) on an unstructured Voronoi mesh grid (see Vogelsberger et al. 2020a, for a review). For self-consistent treatment of ionizing radiation, the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc project employs the radiative transfer (RT) extension A/r.pc/e.pc/p.pc/o.pc-/r.pc/t.pc (Kannan et al. 2019), which solves the first two moments of the RT equation assuming the M1 closure relation (Levermore 1984). The simulation includes the sourcing (from stars and AGNs) and propagation of ionizing photons (in three energy bins relevant for hydrogen and helium photoionization between energy intervals of [ 13 . 6 , 24 . 6 , 54 . 4 , ∞) eV) as well as a non-equilibrium thermochemistry solver to model the coupling of radiation fields to gas. The luminosity and spectral energy density of stars in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc as a complex function of age and metallicity are calculated using the Binary Population and Spectral Synthesis models (BPASS v2.2.1; Eldridge et al. 2017). The sub-grid escape fraction of stars was set to be 𝑓 esc = 0 . 37 in the simulations. This parameter mimics the absorption of Lyman continuum photons happening below the grid scale of the simulation. This parameter was tuned such that the simulated reionization histories approximately match the observed neutral fraction evolution in the Universe. For additional details of the simulation methods, we refer to Kannan et al. (2019, 2022a).  \nIn terms of the galaxy formation model, the simulations employ the IllustrisTNG model (Weinberger et al. 2017; Pillepich et al. 2018a), whichis an update of the Illustris model (Vogelsberger et al. 2014a,b). The simulations include (1) density-, temperature-, metallicityand redshift-dependent cooling of metal-enriched gas (Smith et al. 2008; Wiersma et al. 2009), (2) a two-phase, effective equation of state model for the interstellar medium (ISM) at the sub-resolution level (Springel & Hernquist 2003), (3) star-formation in dense gas following the empirically defined Kennicutt-Schmidt relation, (4) thermal and mechanical feedback from supernovae and stellar winds, (5) metal enrichment from stellar evolution and supernovae, and (6) SMBH formation, growth, and feedback in the quasar and radio modes as described in Weinberger et al. (2017). The model has been extensively tested in large-scale simulations and can produce realistic galaxies that match a wide range of observations, in terms of large-scale galaxy clustering (e.g. Springel et al. 2018), galaxy prop-  \n( 1 ) Name of the simulation. Simulations labeled with ∗ are used only for the numerical convergence study in Appendix A. ( 2 ) 𝐿 box : Side-length of the periodic simulation box. The unit is comoving Mpc (cMpc). ( 3 ) 𝑁 part : Number of particles (cells) in the simulation. In the initial conditions, there are an equal number of DM particles and gas cells. ( 4 ) 𝑚 dm : Mass of DM particles, which is conserved over time. ( 5 ) 𝑚 b : Mass of gas cells in the initial conditions as a reference for the baryonic mass resolution. The gas cells are (de-)refined so that the gas mass in each cell is within a factor of two of this target gas mass. Stellar particles stochastically generated out of gas cells have initial masses typically comparable to 𝑚 b and are subject to mass loss via stellar evolution (Vogelsberger et al. 2013). ( 6 ) 𝜖 : The comoving gravitational softening length for the DM and stellar particles. This is also the minimum gravitational softening length of gas cells, which are adaptively softened. ( 7 ) 𝜖 𝑧 = 6 phy : The physical gravitational softening length at 𝑧 = 6 for reference. ( 8 ) 𝑟 min cell : The minimum physical size of gas cells at the end of the simulations ( 𝑧 ≃ 5 . 5).  \nTable 1. Simulations of the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc suite. Each column corresponds to the following information:  \nerties (e.g. Nelson et al. 2018; Genel et al. 2018; Pillepich et al. 2018b; Naiman et al. 2018), and specifically at high-redshifts (e.g. Vogelsberger et al. 2020b; Shen et al. 2020b, 2022; Kannan et al. 2023).  \nThe initial conditions of the simulations are generated using the G/a.pc/d.pc/g.pc/e.pc/t.pc4 code (Springel et al. 2021) using the second-order Lagrangian perturbation theory at the initial redshift of 𝑧 i = 49. The simulations employ the Planck Collaboration et al. (2016) cosmology with ℎ = 0 . 6774, Ω 0 = 0 . 3089, and Ω 𝑏 = 0 . 0486, 𝜎 8 = 0 . 8159, 𝑛 s = 0 . 9667. We will follow these cosmological parameter choices throughout this paper. In the initial conditions, the gas perfectly follows the DM distribution and is assumed to have primordial composition with hydrogen and helium mass fractions of 𝑋 = 0 . 76 and 𝑌 = 0 . 24, respectively.  \nThe simulations studied in this paper and associated parameters are summarized in Table 1 and all of them are made publically available (Garaldi et al. 2023). Among them, /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 is the flagship simulation of the suite with the largest number of resolution elements in a ( 95 . 5 cMpc ) 3 volume periodic cubic patch of the universe. /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc (Borrow et al. 2022; Shen et al. 2024b) is a subset of high-resolution small-volume simulations using the same numerical setup and physics inputs as the main /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc suite, aiming to explore the formation and evolution of low-mass galaxies in the early Universe. The mass resolution of the L8-N512 (L4-N512) run is about 6 . 4 (51 . 5) times better than the flagship /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 simulation and allows atomic cooling haloes ( 𝑀 halo ∼ 10 8 M ⊙ ; e.g. Bromm & Yoshida 2011; Wise et al. 2014) to be properly resolved with ≳ 200 DM particles, though many important physics mechanisms for haloes below this mass scale have not been included (e.g. explicit treatments of the formation and destruction of molecular hydrogen, Lyman-Werner radiation, and Population-III stars).', '3.1 Galaxy identification and properties': 'We identify dark matter haloes using the Friends-of-Friends (FoF; Davis et al. 1985; Springel 2005) algorithm and further identify inhabiting subhaloes via the S/u.pc/b.pc/f.pc/i.pc/n.pc/d.pc (first described in Springel et al. 2001) algorithm. Gravitationally bound resolution elements are associated with each subhalo, and by definition, exclude elements bound to satellite galaxies. In this paper, galaxies are subhaloes with stellar masses larger than 50 times the baryonic mass resolution ( 𝑚 b ) of the corresponding simulation, including both central and satellite galaxies. A comparison of the sizes of centrals and satellites is presented in Appendix B and displays no significant offsets. Unless  \notherwise stated, we study both central and satellite galaxies without distinction. The centre position of a galaxy is taken from the particle with the minimum gravitational potential energy. The derived galaxy properties that will be used in this paper are summarized in Table 2.', '3.2 Classification of star-forming and quenched galaxies': 'Following common practice in the observational community, we label our galaxies as star-forming versus quenched by recursively searching for the ridge line of the star-formation main sequence (SFMS, defined as the median sSFR of the star-forming galaxies). Following Pillepich et al. (2019), we use 0 . 2 dex bins of stellar mass and reject quenched galaxies as those whose logarithmic distance from the SFMS ridge line at their respective stellar mass lie in the range Δ log 10 sSFR ≤ -1 . 0. We evaluate the SFMS in each stellar mass bin until the number of galaxies in the bin drops below 20 and, in other bins, the reference locus of the SFMS is linearly extrapolated. Star-forming galaxies are those with Δ log 10 sSFR > -1 . 0. This sSFR-based classification is consistent at the 1 -10 percent level with commonly adopted selections in the UVJ diagram (e.g. Fang et al. 2018; Donnari et al. 2019). Due to the dependence of sSFR on both stellar mass and redshift, this selection of star-forming galaxies will be different from a straight cut in sSFR. Furthermore, as discussed in Appendix A, the lower envelope of the SFR𝑀 ∗ distribution could be subject to numerical effects at the low-mass end. Therefore, we conservatively remove the population of spurious quenched galaxies at the low-mass end from our analysis. In Figure 1, we show the identified SFMS at 𝑧 = 6 -10 for /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 and /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc, respectively. Four massive quenched galaxies are found at 𝑧 = 6 and are highlighted as red circles.', '3.3 Mock galaxy images': 'We generate images and datacubes (full spectra in pixels) for simulated galaxies. These images will be used for the comparison with observational results in Section 6. Following our definition of galaxies in Section 3.1, we generate images only for galaxies with stellar mass larger than 50 𝑚 b , and specifically for /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 at 𝑧 ≃ 6 and 8, we increase the stellar mass limit of galaxies to 175 𝑚 b ∼ 10 8 M ⊙ (for /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1) to save computational cost.  \nFollowing Vogelsberger et al. (2020b); Shen et al. (2020b); Kannan et al. (2022b), we use the Monte Carlo dust radiative transfer code S/k.pc/i.pc/r.pc/t.pc (version 8; Baes et al. 2011; Camps & Baes 2015) to generate images of simulated galaxies. We assigned intrinsic emission to stellar particles in the simulations according to their ages and metallicities using the stellar population synthesis method. All stel-  \nFigure 1. The star-formation main sequence (SFMS) of galaxies in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 and /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc at 𝑧 ≃ 6, 8, and 10. The dashed lines are the median relation derived iteratively as discussed in Section 3.2. The four quenched galaxies at 𝑧 ≃ 6 are highlighted with red circles. The orange solid line shows the observationconstrained SFMS at 𝑧 ≃ 6 in Iyer et al. (2020). The purple data points show the compiled SFR𝑀 ∗ of observed 7 ≲ 𝑧 ≲ 13 galaxies from Tacchella et al. (2022, 2023a); Robertson et al. (2023a); Fujimoto et al. (2023).  \n<!-- image -->  \nTable 2. Definitions of halo and galaxy properties in this paper. Properties marked purple are taken directly from the S/u.pc/b.pc/f.pc/i.pc/n.pc/d.pc catalog while others are derived in further post-processing as described.  \nlar particles within a 20 pkpc (physical kpc) aperture are included in this calculation. We adopt the Flexible Stellar Population Synthesis (F/s.pc/p.pc/s.pc) model (Conroy et al. 2009) for continuum stellar emission 1 . For the nebular emission originating from the birth cloud of young star clusters, we use the Byler et al. (2017) model built in F/s.pc/p.pc/s.pc, where line luminosities are generated with C/l.pc/o.pc/u.pc/d.pc/y.pc and tabulated for grids of age, metallicity, gas density, and ionization parameter of the birth cloud. Following Byler et al. (2017), we assume an ionization parameter 2 log 10 𝑈 = -2, which agrees well with the ionization properties of recently JWST -identified high-redshift galaxies (e.g. Reddy et al. 2023). The gas-phase metallicity of the HII regions is chosen to be the same as the initial metallicity of the stellar particles, which is inherited from the gas cell from which a stellar particle is created. Photon packages are randomly released based on the smoothed source distribution characterized by the positions and SEDs of stellar particles. Each stellar particle has been smoothed with smoothing length determined by the distance to the 16th nearest neighbour. Radiative transfer is conducted on a wavelength grid with 657 points spanning from 0 . 05 to 2 𝜇 m, refined around emission lines, as designed in Kannan et al. (2022b).  \nThe emitted photon packages will further interact with the dust in the ISM. We assume that dust is traced by metals in the cold, star-forming gas and turn the metal mass distribution into the dust mass distribution with a constant, averaged dust-to-metal ratio of all galaxies at a fixed redshift. This dust-to-metal ratio depends on redshift as 0 . 9 × ( 𝑧 / 2 ) -1 . 92 , which has been calibrated based on the observed UV luminosity functions at 𝑧 = 2 -10 (Vogelsberger et al. 2020b) and resulted in consistent predictions when applied to /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc (Kannan et al. 2022b). The radiative transfer calculations are performed on an adaptively refined octree grid (Saftly et al. 2014) with the minimum and maximum refinement levels set to 3 and 10. The grid spans 20 pkpc on each dimension and the resulting minimum Octree cell size is about 40 pc. Dust emission is modelled assuming thermal equilibrium with the local radiation field. We adopt the Draine & Li (2007) dust mix, which includes a composition of graphite, silicate, and polycyclic aromatic hydrocarbon (PAH) grains. This model reproduces the average Milky Way extinction curve and is widely used (e.g. Jonsson et al. 2010; Rémy-Ruyer et al. 2014).  \nFor each wavelength on the wavelength grid, 𝑁 p photon packets are launched isotropically. Compared to Kannan et al. (2022b), we increase 𝑁 p to 3 × 10 5 to reduce the Poisson noises in mock images. The photon packets then propagate through the resolved gas (dust) distribution in the ISM and interact with the dust cells randomly before they are finally collected by a pre-defined photon detector, which is placed at a distance of 10 Mpc from the galaxy centre in the positive 𝑧 direction. The pixel-wise galaxy spectra are then recorded, which can be processed to integrated SED, imaging, or integrated field unit data. For rest-frame broadband photometry, galaxy SEDs are convolved with the transmission curves using the S/e.pc/d.pc/p.pc/y.pc code. For the calculation of apparent band magnitudes, the rest-frame flux  \n<!-- image -->  \nFigure 2. Top: Galaxy intrinsic size ( 𝑟 ∗ 1 / 2 ) versus stellar mass at 𝑧 ≃ 6 in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 (blue cloud) and /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc (cyan circles representing individual galaxies). The four quenched galaxies identified in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 are marked with red circles. The vertical dashed line indicates the stellar mass limit we choose for /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. The blue arrow shows the physical gravitational softening length of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. The gray arrows show the convergence radius ( 𝑟 conv ) and the spurious heating mass limit ( 𝑀 lim spur ) of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 (see Appendix A for details). In the massive end, we find that 𝑟 ∗ 1 / 2 anti-correlates with galaxy stellar mass at 𝑀 ∗ ≳ 10 8 M ⊙ . Four compact quenched galaxies are identified in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. Bottom: Redshift evolution of galaxy size-mass relation from 𝑧 ≃ 12 to 𝑧 ≃ 6. The error bars show the median and 1𝜎 dispersion of 𝑟 ∗ 1 / 2 in each bin. The gray shaded region shows the regime closed by 𝑟 conv and 𝑀 lim spur , where numerical effects may prevent the formation of compact galaxies. Negative slopes of the intrinsic size-mass relation are consistently found at these redshifts. Through cosmic time, there is limited evolution at the low-mass end of the size-mass relation while the most massive galaxies continuously become more compact.  \n<!-- image -->  \nis redshifted, corrected for intergalactic medium absorption (Madau 1995), and converted to the observed spectra. The final images have a field-of-view of 15 pkpc and 512 × 512 pixels. Convolution with the point spread function (PSF) of observational instruments or background noises has not been considered here. The pixel size is ∼ 30 pc (physical) and is several times smaller than the observational limit of galaxy size ( ∼ 120 -160 pc) estimated from the half of the FWHM of NIRCam in Morishita et al. (2024).', '4 GALAXY INTRINSIC SIZE-MASS RELATION': 'The intrinsic size of a galaxy is represented by the 3D stellar halfmass radius 𝑟 ∗ 1 / 2 (as defined in Table 2). This is the preferred quantity for our analysis of the physical evolution of galaxy sizes, as it is less affected by geometrical/projection effects and tracer biases (e.g. using light in certain bands as a tracer for stellar mass). Therefore, the analysis will begin with 𝑟 ∗ 1 / 2 and we leave the discussions on observed galaxy sizes to Section 6.1.  \nIn the top panel of Figure 2, we show the intrinsic size-mass relation of galaxies in the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations. The intrinsic sizes of galaxies show a mild increase with stellar mass at the low-mass end ( 𝑀 ∗ ≲ 10 8 M ⊙ ) while displaying a negative correlation with stellar mass towards the massive end. Notably, all galaxies at the massive end ( 𝑀 ∗ ≳ 10 10 M ⊙ ) are compact ( 𝑟 ∗ 1 / 2 ≲ 1 pkpc). Four massive compact galaxies in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 are quenched at 𝑧 ≃ 6. The scatter of galaxy sizes at the low-mass end is substantial. A population of low-mass compact galaxies exists at the same time with extended galaxies with sizes as large as ≳ 5 pkpc. Due to the limited mass and spatial resolution of the simulations, numerical effects can affect galaxy sizes in non-trivial ways (e.g. Power et al. 2003; Ludlow et al. 2019a, 2023). The choice of gravitational softening lengths of the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations is motivated to minimize these effects. This is discussed in detail in Appendix A. In Figure 2, we show two relevant spatial scales for numerical effects, the convergence radius ( 𝑟 conv ) and the physical gravitational softening length ( 𝜖 ) of stars as well as the mass limit below which galaxy size can be affected by spurious heating of DM particles. In this regime, we could miss a population of low-mass compact galaxies in our simulations. However, in terms of the median relation, the numerical convergence is decent when comparing /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 and /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc results. We refer readers to Appendix A for a more thorough discussion.  \nIn the bottom panel of Figure 2, we show the evolution of the intrinsic size-mass relation from 𝑧 ≃ 12 to 𝑧 ≃ 6. The shaded regions show the regime affected by spurious heating at 𝑧 ≃ 6 (as characterized by 𝑟 conv and the spurious heating mass limit). At all redshifts, galaxies at the massive end have their sizes negatively correlated with stellar mass. The massive compact tail continuously evolves towards the massive compact end. As we will discuss later, the tail is represented by galaxies that consistently stay in the compact end of the galaxy population and some of them are quenched at 𝑧 ≃ 6. However, at the low-mass end, there is limited evolution. The scatters of the size-mass relation increase with cosmic time. In Figure 3, we present a more thorough view of galaxy size as a function of stellar mass and redshift. The upper edge of this colormap can be interpreted as the main progenitor mass growth history of the most massive galaxies in the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 simulation. These massive galaxies become increasingly more compact through time since 𝑧 ∼ 10. The compaction phase is similar for slightly lower-mass galaxies but with a delayed starting time. At fixed stellar mass and at 𝑀 ∗ ≳ 10 8 . 5 M ⊙ , galaxy sizes are larger at lower redshift but this is primarily due to their different phases on the track of compaction. On the other hand, the sizes of low-mass galaxies have little evolution and are consistently above ∼ 1 pkpc, suggesting fundamentally different physical mechanisms controlling galaxy sizes at different mass scales.  \nThe negative slope of the size-mass relation at the massive end is consistent with findings in previous theoretical studies using simulations (e.g. Roper et al. 2022; Popping et al. 2022; Costantin et al. 2022; Marshall et al. 2022). However, it is contradictory to the positive correlations found in observations at low redshifts (though at sightly larger stellar masses; e.g. Mosleh et al. 2012; van der Wel et al. 2014; Lange et al. 2015) and some more recent constraints at  \nFigure 3. The evolution median galaxy size as a function of stellar mass and redshift in the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 simulation. The sizes of most massive galaxies become increasingly more compact through time starting at 𝑧 ∼ 10. Slightly lower-mass galaxies follow a similar track but with compaction starting at lower redshifts. On the other hand, the low-mass galaxies have little size evolution with 𝑟 ∗ 1 / 2 consistently above ∼ 1 pkpc.  \n<!-- image -->  \n𝑧 ≳ 6 (e.g. Shibuya et al. 2015; Kawamata et al. 2018; Ormerod et al. 2023; Morishita et al. 2024). The contradiction could be due to the concentrated dust attenuation in massive galaxies, resulting in less cuspy surface brightness profiles and larger apparent sizes (e.g. Roper et al. 2022; Marshall et al. 2022; Popping et al. 2022). In addition, the intrinsic 3D sizes we find in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc are generally larger than the effective radii reported for observed galaxies. A combination of projection effects, tracer biases, and dust attenuation could lead to this discrepancy. This will be discussed in detail in Section 6.2.', '4.1 Comparison to other simulations': 'In Figure 4, we compare the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc results with other cosmological hydrodynamic simulations using diverse numerical setups and galaxy formation physics models. We include predictions from the FLARESsimulation (Vĳayan et al. 2021; Lovell et al. 2021) at 𝑧 = 5, the IllustrisTNG simulation (specifically TNG100-1; Pillepich et al. 2018a) at 𝑧 = 5, the SIMBA simulation (m100n1024; Davé et al. 2019) at 𝑧 = 5, as shown and compared in Roper et al. (2023). In addition, we include predictions of the 2D projected half-stellar-mass radius from the FIRE-2 simulations (Ma et al. 2018b,a) and the Bluetides simulations (Marshall et al. 2022). The 2D half-mass radius is converted to the 3D half-mass radius using a constant multiplication factor, ∼ 1 . 3, as will be motivated in Section 6.2. We also include the UV effective radius predictions of central galaxies from the SERRA simulations (Pallottini et al. 2022), for which we adopt the same multiplication factor and effectively assume that the UV light is a perfect tracer of stellar mass.  \nThe /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc predictions are in good agreement with IllustrisTNG as expected due to the same galaxy formation model employed. Compared to FLARES, /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 predicts roughly the same normalization and negative dependence of size versus mass, but /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 does not predict the sharp galaxy size transition at 𝑀 ∗ ∼ 10 9 . 5 M ⊙ as seen in FLARES. Compared to Bluetides, /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 results are in good agreement with it at the massive end while giving about a factor of  \n<!-- image -->  \nFigure 4. Galaxy intrinsic size-mass relation from a compilation of cosmological hydrodynamic simulations. The /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 and /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc predictions are shown in black open and solid circles, respectively. For comparison, we show results from the FIRE-2 simulations at 𝑧 = 6 , 8 (Ma et al. 2018a), the FLARES simulation at 𝑧 = 5 (Roper et al. 2022, 2023), the IllustrisTNG simulation (Pillepich et al. 2018a) at 𝑧 = 5, the SIMBA simulation (Davé et al. 2019) at 𝑧 = 6, the Bluetides simulations at 𝑧 = 7 (Marshall et al. 2022), and the SERRA simulations (Pallottini et al. 2022). The /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc results are in general consistent with these simulations in shared dynamical range, despite rather different galaxy formation physics models and numerical resolutions of them. An exception is the SERRA simulations which produce increasingly more compact galaxies at 𝑀 ∗ ≲ 10 9 M ⊙ .  \n<!-- image -->  \ntwo larger sizes at 𝑀 ∗ ≲ 10 9 M ⊙ at 𝑧 ∼ 6. Predictions from SIMBA exhibit a positive slope of size-mass relation, in contrast with the results of all other simulations. This could be due to the poorer mass resolution of the m100n1024 run of SIMBA (Roper et al. 2023) but also details in the feedback model.  \nCompared to FIRE-2, which is a suite of high-resolution ( 𝜖 ∼ pc) zoom-in simulations, the results from /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 and /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc are consistent with it at 𝑧 ∼ 6. This consistency is non-trivial given the drastically different simulation strategy, numerical resolution, starformation, and feedback model employed. However, FIRE-2 did not predict the relatively extended galaxy population ( 𝑟 ∗ 1 / 2 ≳ 2 pkpc) at the low-mass end and this discrepancy becomes more apparent at 𝑧 ∼ 8. For IllustrisTNG, the effective ISM model does not re-  \nar phase, and the spatial and time clustering of star-formation may be suppressed compared to simulations employing a high-density threshold and more explicit criteria for starformation (e.g. Pontzen & Governato 2014). Therefore, the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations may miss the temporary extremely compact phase of galaxies, especially at the low-mass end where the stellar component of a galaxy can be dominated by a single compact star cluster (e.g. Yajima et al. 2017; Ma et al. 2018a).  \nEven within high-resolution zoom-in simulations employing resolved ISM physics, predicted sizes can differ substantially due to the detailed implementation of star-formation and feedback. An example is the outlier in our size comparison in Figure 4, the SERRA simulations (Pallottini et al. 2022), which predicts increasingly more compact galaxies at 𝑀 ∗ ≲ 10 9 M ⊙ than all other simulations. As will be discussed later in Section 6, this result turns out to agree better with the observational constraints. These comparisons here highlight the potential of galaxy sizes at high redshifts as a channel to constrain physics models of galaxy formation and evolution.  \nWe note that the simulations compared here together indicate a hump feature in the intrinsic size mass relation of galaxies at 𝑧 ≳ 6 in the stellar mass range of 10 8 to 10 9 M ⊙ . In the mass range, the majority of the galaxies have 𝑟 ∗ 1 / 2 ≳ 1 pkpc. This means that a single power-law relation is not adequate to describe the connection between 𝑟 ∗ 1 / 2 and 𝑀 ∗ at 𝑧 ≳ 6, in contrast to the common practice in observational studies (e.g. van der Wel et al. 2014; Shibuya et al. 2015). Meanwhile, as will be demonstrated in Section 6.1, the extended galaxy sizes around the hump feature are in tension with the latest observational constraints, including the ones based on JWST spectroscopically confirmed galaxy samples.', '4.2 Correlation of size with galaxy properties': 'To understand the physical driver of galaxy size evolution found above, we will first check the correlation of 𝑟 ∗ 1 / 2 with different galaxy properties. In Figure 5, we show the intrinsic size-mass relation of galaxies at 𝑧 ≃ 6 in the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations. The sizes of the square markers are proportional to the logarithm of galaxy number densities in the size-mass plane. The colors map various galaxy properties, including (1) the relative positions of galaxies to the SFMS (top left), (2) the accretion rates of SMBHs normalized by galaxy SFR (top right), (3) the halo gas fraction normalized by the universal baryon fraction ( 𝑓 b ≡ Ω b / Ω m ; middle left), (4) the halo spin (middle right), (5) the Lyman continuum escape fraction ( 𝑓 esc ; bottom left), and (6) the gas disk fraction (bottom right). We remind the readers about definitions of relevant galaxy properties provided in Table 2. In addition to them, 𝑓 esc is defined as the ratio between the number of ionizing photons that escape and are intrinsically emitted from all sources within the virial radius. This is calculated using the method outlined in Yeh et al. (2023). We will focus on at fixed stellar mass , how these secondary parameters affect galaxy sizes. We note that typically in the massive end, the results are affected by small statistics of massive galaxies in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 and should therefore be interpreted with caution.  \nIn the top left panel, we find that galaxies at fixed stellar mass below the SFMS tend to be more compact than the ones above the SFMS. The extreme version of this is the compact quenched galaxies at the massive end, which appear smoothly connected to the compact star-forming population. Similar trends are found at redshift up to ≃ 10. This indicates a potential causal correlation between compaction and quenching of galaxies. This is supported by the trends in the middle left panel, where compact galaxies show lower  \nFigure 5. Galaxy intrinsic size-mass relation at 𝑧 ≃ 6 color-coded by the distance to the SFMS (top left), SMBH accretion rate (top right), gas fraction (middle left), halo spin (middle right), Lyman continuum escape fraction (bottom left), and gas disk fraction (bottom right). The galaxies from /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 are shown by squares and the square sizes are proportional to the logarithm of galaxy density in the size-mass plane. The /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc galaxies are directly shown with circles. The vertical dashed line indicates the stellar mass limit we choose for /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. The shade regions indicate where numerical effects could prevent the formation of most compact galaxies (see Figure 2). Galaxies below the main sequence tend to be more compact and the extreme cases are the compact quenched galaxies at the massive end. At fixed stellar mass, more compact galaxies tend to have more intense SMBH accretion, lower gas abundance, lower halo spin, higher ionizing photon escape fractions, and higher gas disk fraction.  \n<!-- image -->  \n<!-- image -->  \nFigure 6. Galaxy intrinsic size-mass relation at 𝑧 ≃ 5 . 5 color-coded by the matter overdensity (top), and redshift of reionization (bottom) in the nearby environment of galaxies. The notations and plotting styles are the same as Figure 5. To the first order, galaxies with larger stellar masses tend to reside in overdense regions, and their local environments are reionized earlier. Nevertheless, at fixed stellar mass, galaxies in regions that are underdense or have earlier reionization tend to be more compact. However, a population of compact low-mass galaxies resides in extremely dense regions, and their sizes could be significantly affected by environmental effects.  \n<!-- image -->  \nhalo gas mass fractions. In one way, the compact sizes of galaxies could be driven by concentrated starbursts from rapid gas inflow and depletion. Meanwhile, the stellar and AGN feedback associated with the central starburst can drive gas outflows that further decrease the gas abundance and eventually lead to the quenching of galaxies. The causal connection of these signatures will be explored in Section 5 below.  \nIn the top right panel, we find that compact galaxies at 𝑀 ∗ ≳ 10 8 M ⊙ tend to have larger SMBH accretion rates with respect to the galaxy SFR. The SMBH accretion and the concentrated starburst are likely sustained by the same gas inflow towards the centre of the galaxies. We will show a clear relationship between compaction and the ignition of AGN activity in Section 5.1 below.  \nIn the middle right panel, we find that galaxy sizes show a positive correlation with halo spin, which is expected from the classical disk formation theory (e.g. Mo et al. 1998). The role of halo spin in  \ndetermining galaxy size will be investigated in detail in Section 4.4 below. As supporting evidence, in the bottom right panel, we find massive and compact galaxies display prominent disk morphology with more than half of their gas content in a co-rotating component (in terms of angular momentum alignment and within the 2 𝑟 ∗ 1 / 2 aperture as defined in Table 2). More massive and more compact galaxies tend to have higher gas disk fractions. A similar trend exists also for the disk fraction of the stellar component. These findings are qualitatively consistent with the picture of disk instability-driven compaction, as will be discussed in Section 5.3.  \nIn the bottom left panel, we find that the most compact galaxies tend to have orders of magnitude higher Lyman continuum escape fractions than extended galaxies at the same stellar mass. This is consistent with earlier findings in simulations that the gas outflow driven by concentrated starbursts can leave young massive stars (temporarily) transparent to neutral gas absorption and act as the main channel of ionizing photon escape (e.g. Smith et al. 2019; Ma et al. 2020; Barrow et al. 2020; Rosdahl et al. 2022).', '4.3 Correlation of size with environments': 'With the on-the-fly radiative transfer method, the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations can self-consistently track the ionization states of the intergalactic medium around galaxies and enable detailed studies of the relationship between galaxy formation and large-scale environments (e.g. Borrow et al. 2022; Garaldi et al. 2022). We utilize the so-called Cartesian outputs of the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations (Kannan et al. 2022a; Garaldi et al. 2023) at 𝑧 ≃ 5 . 5, which save a subset of gas/dark matter properties on a 3D Cartesian grid. This includes the total matter overdensity ( 𝛿 ) and the redshift of reionization ( 𝑧 reion ), calculated following the method described in Garaldi et al. (2022); Neyer et al. (2023) and Zhao et al. 2024 in prep.. For each halo, we define 𝑧 reion as the latest time at which the (mass-weighted) ionized fraction of the gas in the closest grid cell exceeds the threshold value 0 . 5. For matter overdensity, we take the grid resolution of 512 3 (cell size ∼ 200 ckpc). For 𝑧 reion , we take the grid resolution of 1024 3 (cell size ∼ 100 ckpc). The results have no significant dependence on the choice of grid resolution. For the results presented below, we adopt a smoothing length of 250 ckpc for relevant fields. Results smoothed at protocluster scale (1 cMpc) are presented in Appendix C1 for comparison.  \nIn Figure 6, we show the size-mass relation at 𝑧 ≃ 5 . 5 color-coded by variables related to the environment in which galaxies reside. In the top panel, we show the impact of local matter overdensity. There is a clear leading-order dependence of 𝛿 on galaxy stellar mass. More massive galaxies hosted by more massive haloes tend to reside in denser regions. Nevertheless, at fixed stellar mass , more compact galaxies are more likely to be found in underdense regions while diffuse galaxies preferentially stay in overdense regions. Mergers and gas accretion (associated with denser environments) could play an important role in building extended stellar content of galaxies. This is supported by the larger gas fractions and perturbed morphology of these galaxies as seen in Figure 5. On the other hand, internal physical processes are likely responsible for the formation of compact galaxies. The exception is the population of extremely compact galaxies at 𝑀 ∗ ≳ 10 7 . 5 -10 8 . 5 M ⊙ , which resides in extremely overdense regions. As shown in Figure B1, these galaxies are mostly satellite galaxies and are likely affected by the tidal fields and ram pressure of their hosts. We also note that the impact of matter overdensity shown here is insensitive to the choice of smoothing length. Similar results  \nare found when using a larger 1cMpc smoothing length as shown in Appendix C1.  \nIn the bottom panel, we show the dependence on the redshift of reionization of gas in the environment. Smoothed over the fiducial 250 ckpc scales, there is a strong leading-order dependence that gas around more massive galaxies is reionized earlier in cosmic time. On top of this, at fixed stellar mass , the gas around more compact galaxies tends to experience earlier reionization. The physical interpretation is that, at 250 ckpc scales, local radiation sources and gas budget start to be more relevant for reionization. Massive galaxies that produce more ionizing photons reionize their nearby environment earlier. And at fixed stellar mass , extended galaxies tend to be richer in gas (as seen in Figure 5) and it will take longer time to reionize their local environment given roughly the same amount of ionizing photon budget. We note that, unlike the case of matter overdensity, the dependence of 𝑧 reion does change when using a larger smoothing length as shown in Appendix C. At cMpc scales, 𝑧 reion exhibits a tight positive correlation with matter overdensity. Denser regions reionize earlier due to larger amount of ionizing radiation sources (Garaldi et al. 2022). Therefore, the 𝑧 reion in the size-mass plane closely follows the trend of overdensity as shown in Appendix C1. As a consequence, at fixed stellar mass, 𝑧 reion at cMpc scales is lower around more compact galaxies since they reside in underdense environments. The trend is opposite to what has been found for 𝑧 reion on 250 ckpc scales.  \nAs a summary of the size correlations, in Table 3, we present the Pearson and Spearman rank-order correlation coefficient between 𝑟 ∗ 1 / 2 and the galaxy or environmental properties discussed above.  \nWe perform the measurements around two stellar masses 10 8 and 10 9 M ⊙ . In both mass bins, at least 300 galaxies are included for measuring the correlation coefficients. Out of all the variables considered, the gas fraction has the strongest correlation with galaxy size, followed by halo spin, sSFR, and gas disk fraction. The overdensity shows a strong correlation only at the massive end because of the tip of compact satellite galaxies at the low-mass end. The correlation for overdensity is not sensitive to the choice of the smoothing length. However, for 𝑧 reion , opposite correlations are found when using a smoothing length of 250 ckpc versus 1 cMpc. As discussed above, this reflects a transition from large-scale reionization, which is more affected by overdensity and the amount of external radiation sources, to small-scale reionization, which is more sensitive to local sources and gas abundance. These correlations suggest that (1) a causal connection could exist between compaction, gas depletion, and quenching, (2) external perturbations are likely not the driver for compact but instead generally make galaxies puffier, (3) internal processes related to the disk morphology are likely the trigger of compaction. These aspects will be discussed in more detail in Section 5.', '4.4.1 Disk formation theory': "The correlation of galaxy size with halo spin we see in Figure 5 is expected from the classical theory of disk formation (e.g. Mo et al. 1998; Kravtsov 2013). In this theory, the half-mass radius of the stellar disk is predicted to be  \n𝑟 ∗ 1 / 2 𝑅 vir = 1 . 678 𝑅 d 𝑅 vir = 1 . 678 √ 2 GLYPH<18> 𝑗 d 𝑚 d 𝜆 GLYPH<19> 𝑓 R ( 𝜆, 𝑗 d , 𝑚 d , 𝑐 ) √︁ 𝑓 c ( 𝑐 ) , (1)  \nwhere 𝑅 d is the scale length of the exponential disk, 𝑗 d ( 𝑚 d ) is the angular momentum (mass) ratio of the stellar disk with respect to  \n<!-- image -->  \n<!-- image -->  \nFigure 7. Galaxy intrinsic size versus halo spin at 𝑧 = 6 and 8 in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc1. In each panel, the simulation results are shown by the blue contour with circles and error bars indicating the median relation and 1𝜎 scatter. The colored bands are predictions of an analytical disk formation theory, assuming 𝑚 d = 0 . 01 -0 . 05and 𝑗 d / 𝑚 d values as labelled. The vertical green band shows the prediction from a spherical shell model. In the top two panels, we show the results in mass bin 𝑀 ∗ ≥ 10 9 M ⊙ and 10 8 M ⊙ ≤ 𝑀 ∗ ≤ 10 9 M ⊙ at 𝑧 = 6. In the bottom panel, we show the results in the mass bin 10 8 M ⊙ ≤ 𝑀 ∗ ≤ 10 9 M ⊙ at 𝑧 = 8 (the number of galaxies in the high-mass bin at this redshift is not sufficient enough to form reliable results).  \n<!-- image -->  \nTable 3. Summary of the correlation of intrinsic size ( 𝑟 ∗ 1 / 2 ) with various galaxy/environment properties (in log-log space). The Pearson ( 𝛼 p ) and Spearman rank-order ( 𝛼 s ) correlation coefficients are evaluated based on galaxies with stellar mass log 10 ( 𝑀 ∗/ M ⊙) = 9 ± 0 . 3 or log 10 ( 𝑀 ∗/ M ⊙) = 8 ± 0 . 3. We sort the properties based on the amplitudes of the Spearman correlation coefficients at 𝑀 ∗ ∼ 10 9 M ⊙ . For overdensity and redshift of reionization, we show results when the smoothing length of fields is chosen to be 250 ckpc (fiducial) and 1 cMpc, respectively. Overall, the strongest correlation shows for the gas abundance, followed by halo spin, local matter overdensity, and sSFR with respect to the main sequence.  \nthe host halo, 𝑓 c and 𝑓 R are functions defined in Mo et al. (1998), and the virial radius 𝑅 vir of a halo has been defined in Table 2. We also note the difference between the two definitions of spin ( 𝜆 and 𝜆 ' ) as discussed in Table 2. If we assumed that the baryons and dark matter initially share the same distribution of specific angular momentum (Fall & Efstathiou 1980) and the disk formation process conserves the angular momentum (Mestel 1963) (i.e. 𝑗 d / 𝑚 d ∼ 1), the disk size will be mainly determined by the DM halo properties. In practice, the disk mass fraction 𝑚 d , which enters 𝑓 R , has little impact on the disk size. The halo concentration also has a mild influence on galaxy size. According to the canonical halo mass-concentration relations (e.g. Wechsler et al. 2002; Dutton & Macciò 2014; Diemer & Joyce 2019), high-redshift haloes have 𝑐 ≃ 4 at 𝑧 ≳ 6 with a weak dependence on halo mass. On the other hand, the halo spin parameter 𝜆 was thought to be the dominant factor for disk sizes (e.g. Somerville et al. 2008; Guo et al. 2011; Benson 2012; Somerville et al. 2018).  \nTo quantitatively compare galaxy sizes in the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations with the disk formation theory predictions, in Figure 7, we show 𝑟 ∗ 1 / 2 versus halo spin 𝜆 ' (as defined in Table 2) in different stellar mass bins in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. First, the spin parameter distribution from the simulation is consistent with previous studies using N-body simulations that 𝜆 ' follows a log-normal distribution with median value ∼ 0 . 035 and 1𝜎 scatter ∼ 0 . 5 dex (e.g. Bullock et al. 2001; van den Bosch et al. 2002; Macciò et al. 2007). We also show the predicted 𝑟 ∗ 1 / 2 from disk formation theory, given the median halo mass of galaxies in the stellar mass bin and assuming 𝑚 d in the range 0 . 01 -0 . 05, 𝑗 d / 𝑚 d = 0 . 5 , 1 or 2.  \nFor galaxies with 𝑀 ∗ ≥ 10 9 M ⊙ , they on the bulk agree with the disk formation model with conserved specific angular momentum ( 𝑗 d / 𝑚 d = 1). Galaxies with lower halo spin have smaller sizes. However, we need to note that, within this mass range at fixed halo spin, more massive galaxies tend to have smaller sizes and agree better with model predictions with lower 𝑗 d / 𝑚 d . This could be due to many non-adiabatic processes that can affect the angular momentum of gas fueling star-formation, which will be discussed in Section 5.  \nThe four compact quenched galaxies are typically hosted by low-spin haloes.", '4.4.2 Spherical shell model': "However, galaxies with lower stellar masses (10 8 M ⊙ ≤ 𝑀 ∗ < 10 9 M ⊙ ) display almost no dependence on 𝜆 ' in the low-spin end and only a mild positive dependence in the high-spin end. In the low-mass, low-spin haloes, the geometrical assumption of a stellar disk may break down in the first place. The morphology of the galaxy can be severely disturbed by more frequent mergers and stronger impact of feedback given the shallower potential well of early low-mass haloes (e.g. El-Badry et al. 2018; Jiang et al. 2019).  \nAn alternative proxy for galaxy size is the radius where feedbackdriven outflow encounters the cold gas inflow from CGM and creates a shock. This will be referred to as the spherical shell model. The balance of ram pressure at the shock radius 𝑅 sh gives (Dekel et al. 2023; Li et al. 2023)  \n𝑅 sh 𝑅 str = √︄ 𝑉 w / 𝑀 out 4 𝑉 in / 𝑀 in , (2)  \nwhere 𝑉 in is the gas inflow velocity approximated as the halo virial velocity 𝑉 vir , / 𝑀 in is the inflow rate estimated to be SFR / SFE (SFE is the galaxy-averaged star-formation efficiency, for which we take the empirical parameterization SFE( 𝑀 halo ) outlined in Shen et al. 2023, 2024a). 𝑅 str is the effective radius of accreting gas streams, which can be estimated as (Mandelker et al. 2018; Dekel et al. 2023)  \n𝑅 str ≃ 1 . 66 𝜆 s GLYPH<18> 1 + 𝑧 10 GLYPH<19> 1 / 2 𝑅 vir , (3)  \nwhere 𝜆 s is the contraction factor between the gas stream and dark matter filament and we vary it between 0 . 02 -0 . 14 (Mandelker et al. 2018). 𝑉 w is the wind velocity, which for the IllustrisTNG galaxy formation model adopted by the simulations, has an empirical dependence on halo mass as (Pillepich et al. 2018a)  \n𝑉 w = MIN 𝑉 min w , 𝜅 w × 110 km s -1 GLYPH<18> 𝑀 halo 10 12 M ⊙ GLYPH<19> 1 / 3 GLYPH<18> 𝐻 0 𝐻 ( 𝑧 ) GLYPH<19> 1 / 3 ! , (4)  \nwhere 𝜅 = 7 . 4 and 𝑉 min w = 350 km s -1 . / 𝑀 out is the outflow rate and is related to SFR via the mass loading factor 𝜂 . In the IllustrisTNG model, at low-metallicities ( 𝑍 ≪ 𝑍 w , ref = 0 . 002),  \n𝜂 ≡ / 𝑀 out / SFR = 2 𝑉 2 w ¯ 𝑒 w 𝑁 SNII 𝐸 SNII ( 1 -𝜏 w ) , (5)  \nwhere ¯ 𝑒 w = 3 . 6 and 𝜏 w = 0 . 1 are constants as defined in Pillepich et al. (2018a), 𝑁 SNII = 0 . 0118 M ⊙ -1 is the number of core-collapse supernovae (SNII) per formed stellar mass, 𝐸 SNII = 10 51 erg is the available energy per SNII.  \nIn Figure 7, we show 𝑅 sh with the green band and the value agrees well with the size plateau reached by low-spin galaxies in the stellar mass bin 10 8 M ⊙ ≤ 𝑀 ∗ < 10 9 M ⊙ at both 𝑧 ≃ 6 and 𝑧 ≃ 8. At fixed halo mass, 𝑅 sh has no dependence on halo spin and only a mild dependence on redshift. This is consistent with the slow redshift evolution of galaxy sizes at the low-mass end we found in Section 4. Meanwhile, Equation 2-5 imply that 𝑅 sh ∝ ( 𝑉 w 𝜂 SFE ) 1 / 2 at fixed halo mass. This indicates the potentially strong influence of the feedback prescription on the sizes of low-mass galaxies. Since the wind prescription in the IllustrisTNG model adopted is tightly correlated with halo mass as well, this reduces to 𝑅 sh ∝ SFE 1 / 2 is  \nFigure 8. From top to bottom, we show the evolution of size, gas relative abundance ( 𝑀 gas / 𝑀 ∗ ), sSFR (in unit of Gyr -1 ), and SMBH accretion rate (in unit of M ⊙ yr -1 ) of galaxies from 𝑧 ≃ 13 to 𝑧 ≃ 6. Galaxies are selected at 𝑧 ≃ 6 with 𝑀 ∗ > 10 9 . 5 M ⊙ and are divided into star-forming galaxies and quenched galaxies (see details in the main text). For each galaxy, we follow the main progenitors and document the merger history. The open and solid circles show the times when minor and major mergers happen for the selected quenched galaxies. The evolution of compact quenched galaxies can be divided into three phases, the initial growth in mass and size driven by (minor) mergers, the rapid compaction phase with gas depletion and concentrated star-formation, and the quenching phase.  \n<!-- image -->  \nqualitatively consistent with our findings that galaxy size correlates with their position relative to the main sequence.", '5.1 Evolution of massive compact galaxies': "In this section, we explore the causal connections between galaxy sizes and various galaxy/environment properties found in previous sections. We will first investigate the most compact galaxies at the massive end and follow their evolution history. To achieve this, we subtract the main progenitor branches from the merger trees of galaxies constructed using the LH/a.pc/l.pc/o.pcT/r.pc/e.pc/e.pc algorithm (described fully in the supplementary information of Springel 2005). We track various galaxy properties and identify major (minor) mergers with a merging halo mass ratio threshold of 0 . 1 (0 . 01). The analysis here focuses on massive well-resolved galaxies ( 𝑀 ∗ > 10 9 . 5 M ⊙ ) with unambiguous merger histories. We classify these galaxies into three categories: ex-  \nFigure 9. Similar to Figure 8, here we show the evolution of star-forming galaxies divided into the compact ( 𝑟 ∗ < 1 pkpc) and extended ( 𝑟 ∗ ≥ 1 pkpc) populations at 𝑧 ≃ 6. The notations are the same as Figure 8. We find the compact group tends to have quieter merger histories early on. They later experience stronger gas depletion and have lower sSFR in the end compared to their extended counterparts.  \n<!-- image -->  \n( 𝑟 ∗ ≥ 1 pkpc) star-forming galaxies, compact ( 𝑟 ∗ < 1 pkpc) star-forming galaxies, and quenched galaxies (all turn out to be compact).  \nIn the top panel of Figure 8, we show the time evolution of the sizes of four quenched galaxies selected at 𝑧 ≃ 6, compared to all the star-forming galaxies at a similar mass scale. In the bottom three panels, we show the time evolution of the galaxy gas abundance (normalized by stellar mass, 𝑀 gas / 𝑀 ∗ ), sSFR, and SMBH accretion rates of the same galaxies, respectively. The major and minor mergers of the quenched galaxies are shown by solid and open circles. We identify three phases of the evolution of a compact quenched galaxy in the first Gyr of the universe. The later quenched population is initially indifferent from the star-forming galaxies with stellar mass and size growth dominated by minor mergers and smooth accretion. During this epoch, the baryon content of the galaxies is dominated by gas and galaxies grow in a self-similar fashion maintaining a stable sSFR. In the second phase beginning around 𝑧 ∼ 10, rapid gas depletion and compaction of the stellar content happen at the same time in these later quenched galaxies. Right after compaction starts, the sSFRs of these galaxies rise temporarily compared to the normal starforming counterparts followed by a mild decline. During this phase, SMBH accretion gets activated at the same time as concentrated star-formation happens. None of these phenomena are associated with major or minor merger events, suggesting internal processes dominating the size evolution. In the final phase, the sSFR drops in quenched galaxies in an extremely short time scale ( ≲ 10 Myr) and so does the SMBH accretion. The gas is completely depleted and the sizes of these quenched galaxies stabilize at ≲ 0 . 5 pkpc later on.  \nIn Figure 9, we show the same type of evolution history for the star-forming galaxies at 𝑧 ≃ 6, divided into compact and extended  \nFigure 10. SFR surface density versus stellar mass density of galaxies in the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 simulation compared to observational samples (Morishita et al. 2024). We highlight the evolution trajectories of compact quenched galaxies from 𝑧 ≃ 12 to 𝑧 ≃ 6. Each marker is separated by Δ 𝑧 = 1. The arrow indicates the critical surface mass density threshold motivated by previous theoretical studies of massive compact galaxies (known as 'blue/red nuggets', BN/RN) at cosmic noon.  \n<!-- image -->  \npopulations. Within the star-forming galaxies, we find a distinction between the compact and the extended populations and the trend is similar to what we found for quenched and star-forming galaxies above. Starting around 𝑧 ∼ 10, the compact star-forming galaxies have comparably lower gas abundance as their sizes deviate from the extended population, followed by declined sSFRs later on. These findings are broadly consistent with the trends we found in Figure 5 that compact galaxies tend to have lower sSFR and gas abundance. These compact star-forming galaxies might have gone through a similar gas depletion and compaction process to the quenched population but to a milder extent or at shifted timings due to lower masses or higher spins. Compared to the extended ones, the compact starforming galaxies also have quieter merger histories, which is consistent with the underdense environments they reside in as found in Figure 6. This again suggests that internal processes are responsible for the compaction of galaxies rather than external perturbations.  \nThe rapid quenching of the few compact galaxies is likely driven by the rapid removal of cold gas reservoirs by stellar and AGN feedback. Evidences are shown in the gas images in Figure 11, where the central gas distribution of the quenched galaxy is disturbed and displays shock-like features. However, we are not able to disentangle the relative contribution of stellar and AGN feedback in quenching. In addition, the development of a compact bulge can help stabilize the disk instability and drive morphological quenching (e.g. Martig et al. 2009). However, this may not completely shut down star-formation in the entire galaxy but only prevents mass inflow and central starburst. Other external mechanisms for quenching will operate on longer time scales than we observe in these simulations. For example, secular halo quenching can happen due to the development of a hot virialized inner CGM, which suppresses cold gas fueling the galaxy. These could be relevant for the mild decrease in sSFR in compact star-forming galaxies.  \nIn Figure 10, we show the SFR surface density ( Σ SFR ≡ SFR /( 2 𝜋 𝑟 ∗ 2 1 / 2 ) ) versus stellar mass surface density ( Σ ∗ ≡  \n𝑀 ∗ /( 2 𝜋 𝑟 ∗ 2 1 / 2 ) ) of galaxies in the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 simulation. They are compared to the observational constraints from JWST (Morishita et al. 2024). We highlight the evolution trajectories of the quenched galaxies from 𝑧 ≃ 12 to 𝑧 ≃ 6. Σ SFR follows a tight correlation with Σ ∗ with little redshift evolution and shows good agreement with the observed sample despite lower scatter. The later quenched galaxies grow their stellar mass following this correlation until reaching a stellar mass surface density threshold and become quenched afterward. This threshold is about Σ ∗ ∼ 10 9 -10 10 M ⊙ pkpc -2 , corresponding to Σ SFR ∼ 10 1 -10 2 M ⊙ yr -1 pkpc -2 . This threshold happens to be consistent with the compaction and quenching phases discussed extensively in literature for galaxies at cosmic noon (e.g. Dekel et al. 2009b; Dekel & Burkert 2014; Zolotov et al. 2015; Tacchella et al. 2016; Lapiner et al. 2023) as well as the early proto-bulge formation found in observations (e.g. Baggen et al. 2023; Baker et al. 2023).", '5.2 Gas angular momentum': "To understand the physical reason for compaction, we need to start from the beginning when the gas fueling star-formation is accreted into the halo. The angular momentum of the gas has long been considered one of the most important factors for regulating galaxy sizes. We can gain insights into the time evolution of gas angular momentum by examining the specific angular momentum (or equivalently spin) profile of the halo.  \nIn Figure 12, we show the spin profile of galaxies in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 at 𝑧 ≃ 6 grouped in two stellar mass bins. The spin of a particle species 𝑥 is defined as 𝑗 𝑥 /( √ 2 𝑉 vir 𝑅 vir ) , where 𝑗 𝑥 is the specific angular momentum of the corresponding particle species. This definition is chosen such that, for dark matter, it becomes equivalent to the halo spin defined in Table 2. Here we consider the spin profile of dark matter, stars (restricted to 𝑟 < 2 𝑟 ∗ 1 / 2 ), all gas, and specifically cold gas  \n( 𝑇 < 10 5 K). The spin of dark matter increases nearly linearly with radius and the central angular momentum has little correlation with the outskirts. Similar spin or specific angular momentum profiles have been found in previous cosmological 𝑁 -body simulations (e.g. Bullock et al. 2001; Bett et al. 2010) and is consistent with a toy model for halo buildup by tidal stripping of satellites as they spiral inward via dynamical friction (Bullock et al. 2001; Maller et al. 2002).  \nIn the classical theory, haloes acquire angular momentum from the gravitational tidal torques exerted by neighbouring density perturbations (e.g. Hoyle 1951; Peebles 1969; Doroshkevich 1970; White 1984; Barnes & Efstathiou 1987) and through mergers (e.g. Vitvitska et al. 2002; Maller et al. 2002). Baryons were expected to initially follow the dark matter mass distribution and share the same specific angular momentum distribution. This motivates the classical disk formation theory described in Section 4.4. In this picture, halo gas goes through a joint semi-spherical infall with dark matter from outside the halo before being heated to 𝑇 vir by a virial shock and becoming pressure-supported. The halo gas then cools down, contracts to the bottom of the halo potential while conserving angular momentum, and eventually forms a rotation-supported disk (e.g. Fall & Efstathiou 1980; Mo et al. 1998; Bullock et al. 2001).  \nHowever, at high redshifts, halo gas accretion is dominated by the 'cold mode' (e.g. Kereš et al. 2005; Dekel & Birnboim 2006; Kereš et al. 2009; Dekel et al. 2009a; Wetzel & Nagai 2015) driven by supersonic cold gas streams penetrating the halo directly to the central disk (often spatially associated with dark matter filaments) as well as gas clumps in merging subhaloes. The cold gas distribution shown in Figure 11 clearly demonstrates these features. Therefore,  \nFigure 11. Surface density maps of the moderately cold ( 𝑇 < 3 × 10 5 K) gas distribution around two simulated galaxies selected at 𝑧 ≃ 6 and their main progenitors out to 𝑧 ≃ 12. The first galaxy is a compact quenched galaxy while the second is an extended star-forming galaxy. The cold gas distributions in the CGM are dominated by filament-like gas streams. The star-forming galaxy has a central peak in the cold gas distribution while the cold gas distribution of the quenched one is disturbed by feedback at the centre.  \n<!-- image -->  \ndeviations of the specific angular momentum of gas and dark matter are expected (e.g. Kereš et al. 2009; Stewart et al. 2013; Nelson et al. 2015). At the outskirt of the halo ( 𝑟 ∼ 𝑅 vir ), the specific angular momentum of initially accreted cold gas can be larger than the dark matter due to stronger quadruple torques experienced by the gas (Danovich et al. 2015), but this difference will be weak at high redshifts. In the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations, we do find a small enhancement of the spin of cold gas compared to dark matter at the halo outskirts. Meanwhile, we also find that the gas spin is generally higher than that of the dark matter component at 𝑟 ≲ 0 . 5 𝑅 vir and often features a plateau-like feature until reaching 𝑟 ∼ 0 . 1 𝑅 vir . This is consistent with many earlier studies (e.g. Sharma & Steinmetz 2005; Kimm et al. 2011; Danovich et al. 2015; Teklu et al. 2015). This is because the newly accreted cold gas has not yet mixed with the earlier accreted ones that have lower angular momentum and cluster in the central region due to cooling. The inflowing gas streams roughly conserve their angular momentum (Stewart et al. 2013; Danovich et al. 2015). The plateau of gas spin happens to be similar to the halo average spin 𝜆 ' even for the relatively compact population. This suggests that there is enough angular momentum and rotation support of accreted gas before reaching the vicinity of the galaxy, and it should be the physical processes within the galactic radius that drive the compaction of the stellar content.  \nIn the messy region below the impact parameter of the streams, ≲ 0 . 1 𝑅 vir , tidal torques from the inner galactic structure as well as complex interactions of streams with other streams and feedbackdriven outflows can substantially decrease the angular momentum of gas (e.g. Ceverino et al. 2010; Danovich et al. 2015; Nelson et al. 2015). Closer to the star-forming region, the spin of gas starts be aligned with the stellar disk. Further loss of gas angular momentum at smaller radii is driven by the disk instability (e.g. Dekel et al. 2009b; Dekel & Burkert 2014, discussed later in Section 5.3). These processes together result in the decline of gas spin profile at small  \ngalactocentric radii, which can be substantially smaller than the overall averaged halo spin. Disk instability is more likely to be the driver of the compaction of massive galaxies where feedback-driven outflows are less effective and disk morphologies show up.  \nAt similar radial ranges, we find the spin alignment of gas and stars becomes poorer below 𝑟 ∗ 1 / 2 , which is a consequence of the proto-bulge formation from concentrated star-formation. We also find that the spin of gas is about 2 -3 times larger than that of the stellar disk. Similar gaps have been found in observational studies as well (Obreschkow & Glazebrook 2014; Cortese et al. 2016). As a more direct tracer of the recently arrived mass, gas is expected to come in with a higher specific angular momentum, while stars are made of baryons that arrived earlier with lower specific angular momentum (e.g. Danovich et al. 2015; Teklu et al. 2015). In addition, feedback from young massive stars and potentially AGN can drive galactic outflows, which preferentially remove the low-angular momentum gas in the dense central regions of the galaxy and effectively increase the specific angular momentum of the remaining gas (e.g. Maller et al. 2002; Governato et al. 2010; Brook et al. 2011, 2012; Guedes et al. 2011; Romanowsky & Fall 2012; Agertz & Kravtsov 2016; Lagos et al. 2017). Such effects are likely more prominent at high redshifts, where star-formation is more clustered and bursty (e.g. Smit et al. 2016; Sparre et al. 2017; Emami et al. 2019; Flores Velázquez et al. 2021; Shen et al. 2023; Sun et al. 2023b).  \nCompared to galaxies in the low-mass bin, massive galaxies have better alignments of gas and stellar angular momentum at 𝑟 ≲ 0 . 1 𝑅 vir and clear signatures of a co-rotating stellar disk at 𝑟 ∗ 1 / 2 ≲ 𝑟 ≲ 2 𝑟 ∗ 1 / 2 . This is consistent with the fact that massive galaxies are better described by the classical disk formation theory while low-mass galaxies are better described by the spherical shell model shown in the section above. This change in geometry could be due to the more frequent mergers and stronger impact of feedback-driven outflows  \nFigure 12. Top : Spin profile of dark matter, gas, and star at 𝑧 = 6 in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. We show the stacked results of galaxies in the mass bin of 𝑀 ∗ ≥ 10 9 M ⊙ (left) and 10 8 M ⊙ ≤ 𝑀 ∗ ≤ 10 9 M ⊙ (right) . We show the results of all galaxies in the mass range in solid and the compact ( 𝑟 ∗ 1 / 2 < 1 pkpc) subset of galaxies in dashed lines. The gray circle with the error bar indicates the median spin of dark matter haloes ⟨ 𝜆 ' halo ⟩ ≃ 0 . 035. The red circle with the error bar indicates the median and 1𝜎 scatter of 𝑟 ∗ 1 / 2 of stacked galaxies. Bottom : Spin alignment of dark matter, gas, and star to the galaxy face-on direction defined based on stellar content (as defined in Table 2). 𝜃 is the angle between the spin vector of each component with the galaxy face-on direction.  \n<!-- image -->  \nof low-mass galaxies. Compared to the overall galaxy population, compact galaxies tend to have lower spins of both dark matter and gas upon accretion but similar spin profiles in the inner region of the halo. The initial lower spin makes these haloes more prone to disk instability (Dekel & Burkert 2014) and enters the compaction phase earlier.", '5.3 Disk instability': 'One important driver for the compaction of massive galaxies at 𝑧 ≳ 6 is disk instability. As shown in Figure 5 and Figure 12, we find more prevalent signatures of disk morphology for both the stellar and gas components of these galaxies. This is expected from disk formation theory as massive galaxies tend to have deeper and more concentrated gravitational potentials (e.g. Hopkins et al. 2023) and increased pressured confinement from the hot virialized CGM (e.g. Stern et al. 2021). However, the gaseous disks in high-redshift starforming galaxies are likely not stable or in a marginally stable state.  \nThe large cold gas abundance and rapid gas cooling can trigger global instability in the disk. This can be understood through a modified turbulent Toomre instability 3 criterion of a differentially rotating', 'disk (Toomre 1964)': "𝑄 = √︃ 𝑐 2 s + 𝜎 2 r 𝜅 𝜋 𝐺 Σ tot > 𝑄 c , (6)  \nwhere 𝜅 is the epicyclic frequency (e.g. √ 2 Ω for a flat rotation curve, where Ω is the angular circular velocity), 𝑐 s is the thermal sound speed, 𝜎 r is the radial velocity dispersion, Σ tot is the total matter surface density, and 𝑄 c is an order-unity constant. Disk instability happens if the LHS drops below 𝑄 c , which can be a consequence of either cooling or increased Σ tot from gas accretion. Note that even when turbulence dominates over the thermal support, turbulence cascade to the viscous scale can still turn turbulent energy into thermal energy, which is eventually carried away by cooling radiation.  \nIn high-redshift galaxies, both scenarios are highly possible. For example, in the standard Λ CDMcosmology, the typical gas accretion time scale can be estimated as 𝑡 acc = 𝑀 halo / / 𝑀 halo ≃ 316 Myr ( 1 + 𝑧 / 7 ) -2 . 25 ( 𝑀 halo / 10 12 M ⊙ ) -0 . 15 (Dekel et al. 2009b), where the gas accretion rate of a halo has been assumed to scale with the total mass growth rate by a constant factor. The ratio of 𝑡 acc to the Hubble time 𝑡 H ≡ 1 / 𝐻 ( 𝑧 ) will become lower at higher redshift and can drop significantly below unity at 𝑧 ≳ 6 to allow instabilities to develop in disks. A similar thing applies to the typical cooling time in the starforming disk (e.g. Dekel et al. 2023) which is almost always below ∼ 100 Myr.  \nFigure 13. Images of two galaxies at 𝑧 ≃ 6 and 𝑧 ≃ 8 in the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 simulation in rest-frame UV and optical (V band) with and without dust attenuation. The images are generated following the procedure described in Section 3.3. The orange circles indicate the effective radius in the corresponding band centered on the surface brightness-weighted centre of the galaxy.  \n<!-- image -->  \nOnce instabilities get triggered, clumps with sizes comparable to the Jeans length will form over a dynamical crossing time scale  \n𝑡 d ≡ 1 Ω = 𝑅 d 𝑉 c ≃ 10 Myr GLYPH<18> 𝑅 d 1 pkpc GLYPH<19> GLYPH<18> 𝑉 c 100 km s -1 GLYPH<19> -1 , (7)  \nwhere 𝑅 d ∼ 𝑟 ∗ 1 / 2 is the characteristic size of the disk and 𝑉 c is the circular velocity at 𝑅 d . These clumps will migrate toward the centre of the galaxy through clump-clump interactions and meanwhile drive general mass inflows 4 . In a self-regulated state, Dekel et al. (2009b) showed that the evacuation time scale of the disk due to clump migration is  \n𝑡 evac ≃ 20 𝑡 d GLYPH<16> 𝛼 0 . 2 GLYPH<17> -1 GLYPH<18> 𝑄 c 0 . 68 GLYPH<19> 2 GLYPH<18> 𝛿 0 . 5 GLYPH<19> -2 , (8)  \nwhere 𝛿 is the disk mass to total mass ratio within 𝑅 d , 𝛼 is the disk mass locked in giant gas clumps as defined in Dekel et al. (2009b), and 𝑄 c = 0 . 68 for an isothermal thick disk (Goldreich & LyndenBell 1965). The compaction happens faster in gas-rich environments with high 𝛿 . Since baryons dominate the mass budget in the central halo and gas dominates over stars in the early phase of galaxy growth (see Figure 8), 𝛿 is comparable to the 𝑓 gas disk measured in Figure 5 and can reach ≳ 0 . 5 in massive galaxies on the track of compaction. The time scale for instability-driven compaction will be of order a few hundred Myr, which is quite consistent with the evolution track shown in Figure 8. As argued in Dekel & Burkert (2014), once disk instability gets triggered, 𝛿 will monotonically increase during the gas contraction phase and this will be a runaway process. This could  \ncreate a bimodal track of galaxy size evolution that depends on halo spin, which has been discussed extensively in the context of compact galaxies at 𝑧 ∼ 2 -3 (e.g. Dekel & Burkert 2014; Zolotov et al. 2015; Lapiner et al. 2023). However, at higher redshifts, the compaction track will quickly become much more likely to happen (Dekel & Burkert 2014).  \nThe gas inflow is also likely to be 'wet' (Dekel & Burkert 2014; Krumholz et al. 2018). In the self-regulated marginally stable disk, star-formation happens at the time scale 𝑡 sf = 𝑡 ff / SFE ff ≃ 50 𝑡 d ( SFE ff / 0 . 01 ) -1 , where SFE ff is the star-formation efficiency in molecular phase per free-fall times and we take the value ∼ 0 . 01 suggested by observations (e.g. Sun et al. 2023a, needs to be distinguished from the galaxy-scale SFE defined above), 𝑡 ff ≃ 0 . 5 𝑡 d is the free-fall time assuming star-formation happens in regions that are denser than the mean density of the disk by a factor of few (Dekel & Burkert 2014). 𝑡 sf is usually longer than 𝑡 evac , especially when 𝛿 is large, cold gas inflow can be sustained against local star-formation in the disk, and therefore, fuels the starburst at galaxy centres and helps form compact galaxies.", '6.1 Apparent size-mass(luminosity) relations': "The analyses above focus on the intrinsic physical sizes of galaxies. However, in the observations, galaxy sizes are often measured by fitting the surface brightness profile/distribution of galaxies in certain bands and quoted in terms of the effective radius ( 𝑅 eff , 2D 'halflight' radius). To enable a fair comparison of simulation results with observations, we utilize the mock galaxy images described in Section 3.3. We find 𝑅 eff as the circular aperture centered on the  \n<!-- image -->  \n<!-- image -->  \nFigure 14. Apparent size-mass relation at 𝑧 ≃ 6, 8, 10. The open circles and error bars show the median and 1𝜎 scatter of galaxy effective radii ( 𝑅 eff ) in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. Different colors correspond to the size in UV, V band, and UV without dust attenuation. The background colormap shows the full distribution of 𝑅 eff in UV versus stellar mass. The solid circles show individual galaxies in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc on this plane. The vertical dashed line indicates the stellar mass limit we choose for /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 mock images. We compare simulation results with observational constraints from Morishita et al. (2024) as shown by open gray squares with error bars.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 15. Apparent size-luminosity relation at 𝑧 ≃ 6, 8, 10. The open circles and error bars show the median and 1𝜎 scatter of galaxy UV effective radii ( 𝑅 UV eff ) in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 with and without dust attenuation. The background colormap shows the full distribution of 𝑅 UV eff versus 𝑀 UV . The solid circles show individual galaxies in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc on this plane. We compare simulation results with observational constraints from Shibuya et al. (2015); Kawamata et al. (2018); Yang et al. (2022); Ono et al. (2023a,b).  \n<!-- image -->  \npixel flux-averaged centre 5 of the image that contains half of the total light of all image pixels. In Figure 13, we show examples of these mock galaxy images in rest-frame UV and optical V band with and without dust attenuation. The galaxy 𝑅 eff in UV and V bands are fairly close to each other before dust attenuation. However, dust attenuation increases the apparent size of galaxies preferentially in UV, which will be discussed in detail below.  \nIn Figure 14, we show the galaxy 𝑅 eff versus stellar mass relations in rest-frame UV and V band in the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations. We also show the UV 𝑅 eff when dust attenuation is not taken into account. Thediscrepancies of size from band choice or dust attenuation mainly show up in massive galaxies. Unlike the intrinsic 𝑟 ∗ 1 / 2 , the 𝑅 eff in UV display a mild positive correlation with galaxy stellar mass. The simulation results are compared to the measurements in Morishita et al. (2024), which were based on 341 galaxies at 5 < 𝑧 < 14 in a suite of nine public JWST extragalactic fields 6 in Cycle 1. We find that the observed galaxies are systematically more compact than the simulated ones. A large population of compact galaxies exists at 𝑀 ∗ ≲ 10 9 M ⊙ with 𝑅 eff ≲ 200 pc. We note that they are still larger than the full width at half maximum (FWHM) of the PSF of short-wavelength NIRCam filters (e.g. for F070W, F090W, F115W, FWHM ∼ 0 . 03 arcsec ∼ 150 pc at 𝑧 ∼ 8; Rigby et al. 2023).  \nIn Figure 15, we show the galaxy 𝑅 eff -𝑀 UV relation from the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations. The absolute UV magnitude 𝑀 UV is derived by integrating the fluxes of galaxy images in rest-frame UV. Similar to the size-mass relation, the size-luminosity relations after dust attenuation show a mildly positive slope at 𝑧 ≃ 6 and are almost flat at 𝑧 ≳ 8. We compare the simulation results with observational data based on HST from Shibuya et al. (2015) and Kawamata et al. (2018) as well as the ones based on JWST from Yang et al. (2022); Ono et al. (2023a,b). Similar discrepancies between simulations and observations are found at the faint end ( 𝑀 UV ≳ -20) at about 1𝜎 level.  \nDue to the numerical limitation, the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations could miss a population of compact galaxies below the convergence scale and the gravitational softening length (as discussed in Section 4), which is indicated by the gray shaded region in the figure. However, the offset between the median relations is robust against numerical effects when taking the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc results into account. On the other hand, the discrepancy in size could indicate limitations of the feedback model employed by the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations (and also IllustrisTNG). We note this is not restricted to /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc, but essentially many other simulations shown in Figure 4 fail to produce compact ( 𝑅 eff ≲ 300 pc) galaxies at the stellar mass ≲ 10 9 M ⊙ at 6 ≲ 𝑧 ≲ 10.  \nIn Ono et al. (2023b), a population of extremely compact ( 𝑅 eff ≲ 100 pc) galaxies was compared to the FOREVER22 (Yajima et al. 2022) and Renaissance simulations (Xu et al. 2016; Barrow et al. 2017), which produced a few analogs at higher redshifts ( 𝑧 ≳ 10) in overdense regions. These galaxies were argued to be formed under the impact of galaxy mergers or interactions. This scenario is similar  \n5 The median offset between the flux-average centre and the galaxy centre defined using potential minimum in Section 3.1 is ∼ 500 physical pc, and, in ≳ 80% of the cases, the offset is ≲ 1 pkpc. This is not a negligible effect given the small apparent size of some of these galaxies.  \n6 The fields included in Morishita et al. (2024) analysis are PAR1199 (Stiavelli et al. 2023, PID 1199;), J1235 (PID 1063; PI: Sunnquist), NEP (Windhorst et al. 2023, PID 2738;), PRIMER-UDS and COSMOS (PID 1837; PI: Dunlop), NGDEEP (Bagley et al. 2024, PID 2079; PI: Finkelstein;), CEERS (Bagley et al. 2023; Finkelstein et al. 2023, (PID 1345;), GLASSJWST/UNCOVER(Treuetal. 2022, PID 1324;), and JADESGDS (Eisenstein et al. 2023; Robertson et al. 2023b; Tacchella et al. 2023a).  \nto the subpopulation of compact galaxies at 𝑀 ∗ ≲ 10 8 M ⊙ we find in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc in extremely overdense regions (see Figure 6). However, they are not representative of the general galaxy population. In contrast, we find that the sizes of the majority of the low-mass galaxies increase during the phase of merger-driven growth and are positively correlated with matter overdensity.  \nIn the SERRA simulation (Pallottini et al. 2022), compact galaxies were consistently produced at 𝑧 ≃ 7 . 7 with 𝑅 UV eff ∼ 100 pc at 𝑀 ∗ ∼ 10 9 M ⊙ and were on the other hand slightly below the observational constraints. In the high-resolution variants of the S/i.pc/m.pc/b.pc/a.pc simulations (Davé et al. 2019; Wu et al. 2020), galaxy sizes were found to decrease down to ∼ 200 pc at 𝑀 UV ∼ -18 at 𝑧 ≃ 6. These dramatic and qualitative differences between size predictions of different simulations may indicate incompleteness in many current feedback implementations. The prescriptions configured to prevent excessive star-formation at a similar level might result in rather different kinematic behavior of the gas that fuels star-formation (Rosdahl et al. 2017; Hopkins et al. 2018) and different stellar morphology. This highlights the potential of using galaxy morphology to constrain galaxy formation models at high redshifts.", '6.2 Observational factors shaping apparent sizes': 'In addition to the theoretical uncertainties discussed above, several observational effects could be relevant and cause order unity differences between the intrinsic and observed sizes. We will discuss them in this section. In Figure 16, we show the ratios of galaxy sizes with various definitions. The projection of the galaxy is taken to be the positive 𝑧 direction in the simulation coordinates and effectively represents a random viewing angle when analyzing an ensemble of galaxies. We summarize several key effects that set the difference between the observed galaxy size in certain band (usually rest-frame UV) and the intrinsic half-mass radius.  \nProjection effects : Since a spherical aperture always contains less mass or light than a cylindrical aperture, projection tends to make galaxies appear smaller. We find that the ratio between the 2D halfmass radius (with random projection) and the 3D one is insensitive to galaxy stellar mass and takes a median value of ∼ 0 . 78. This roughly agrees with the estimates in Wolf et al. (2010) assuming spherical symmetry and a variety of surface brightness profiles. This also motivates the multiplication factor we used in comparing different simulations in Section 4.1.  \nTracer effects : We compare the 𝑅 eff in UV (before dust attenuation) with the 2D half-mass radius under the same viewing angle. We aim at understanding how well UV light acts as a tracer of the stellar mass in observations. We find that, at 𝑀 ∗ ≲ 10 9 M ⊙ , the UV effective radius is smaller than the stellar mass effective radius and differs by a constant factor of ∼ 0 . 8. This is because recent starbursts take place in the central part of the galaxy while the old and potentially ex-situ stars reside in the outer part. At 𝑀 ∗ ≳ 10 9 M ⊙ , the ratio gradually increases to unity due to the compaction of the galaxies and that most of the stellar mass is dominated by the recent central starburst.  \nDust attenuation : We find that the UV effective radius after dust attenuation is larger than that before dust attenuation. The ratio between the two increases monotonically with stellar mass and can reach ≳ 2 at 𝑀 ∗ ≳ 10 10 M ⊙ . This is mainly due to the relatively large abundance of dust in dense ISM and their spatial correlation with the young stellar populations. Therefore, dust attenuation is more concentrated in the central region of the galaxy and smooth out the surface brightness profiles of galaxies. As seen in Figure 14, it qualitatively changes the slope of the 𝑅 eff -𝑀 ∗ relation from  \nFigure 16. An illustration of key factors in bridging the apparent and intrinsic galaxy size at 𝑧 ≃ 6. From top left to bottom right, we show the size ratios between (1) 𝑅 eff in UV after and before dust attenuation, (2) 𝑅 eff in UV and V band, (3) 𝑅 eff in UV and 𝑅 eff of the stellar mass, (4) 𝑅 eff of stellar mass and the intrinsic 3D 𝑟 ∗ 1 / 2 , respectively. The vertical dashed line indicates the stellar mass limit we choose for /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 mock images. Dust attenuation makes galaxies appear larger in UV and this effect is more prominent in more massive galaxies. The 𝑅 eff in UV is slightly larger than at longer wavelengths mainly due to the difference level of dust attenuation. The 𝑅 eff in UV is slightly smaller than the 𝑅 eff of underlying stellar mass distribution due to more young stars forming in the central part of the galaxy but this difference diminishes at the massive end. The 2D sizes of galaxies are consistently smaller than the 3D sizes by about 20 percent due to projection effects.  \n<!-- image -->  \nintrinsically negative to positive, which is more consistent with observations. Similar effects have been discussed extensively in many earlier works (e.g. Wu et al. 2020; Roper et al. 2022; Marshall et al. 2022; Popping et al. 2022; Cochrane et al. 2023).  \nObserving band : We find a mild ( ∼ 10 -30 percent) increase in galaxy 𝑅 eff in UV compared to the optical V band. The tracer effect discussed above could have a minor impact on the sizes at different wavelengths. However, the trend here is primarily driven by stronger dust attenuation in rest-frame UV. At the mass range around 𝑀 ∗ = 10 8 M ⊙ , there is almost no difference between the two since the dust attenuation is minimal in low-mass galaxies. This is in good agreement with observational results at a similar mass scale (Morishita et al. 2024).  \nSurface brightness bias : Given the background noise, observations could completely miss sources if their peak surface brightness is below the detection limit. We find some signatures of this in Figure 10 as the observed sample cuts off at Σ SFR ∼ 1 M ⊙ yr -1 pkpc -2 . This can be translated to the UV surface brightness (using the same Kennicutt 1998 relation adopted in Morishita et al. 2024) 𝜇 lim ∼ 28 mag arcsec -2 , which is close to the detection limit of  \ncurrent JWST surveys 7 . The impact of this incompleteness has been discussed in Kawamata et al. (2018) and can lead to a ∼ 0 . 2 -0 . 3 dex increase in median galaxy sizes at 𝑧 ∼ 6 -9 as well as increasing the steepness of the size-luminosity relation.  \nEven for identified sources, the surface brightness profile can only be probed out to a certain radius beyond which extrapolations of the Sérsic profile are required. This can introduce a bias as the fits derived on the central region may fail dramatically at the wings and lead to artificially smaller apparent sizes for intrinsically diffuse galaxies (Roper et al. 2022). In the top panel of Figure 17, we show the stacked surface brightness profile of 200 galaxies with 𝑀 ∗ ∼ 10 9 M ⊙ at 𝑧 ≃ 6 in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. After stacking, the surface brightness distribution is roughly axially symmetric and we measure the surface brightness profile in a set of linearly-spacing annuli. We then fit the  \n7 We perform a simple estimate by taking the 5𝜎 point source limiting magnitudes ∼ 29 mag in e.g. Morishita et al. (2024); Ormerod et al. (2023), dividing the corresponding flux limit by the area of the assumed aperture 𝑟 0 = 0 . 16 arcsec, and converting it to 1𝜎 limit. We obtain 𝜇 lim ≃ 29 + 2 . 5 log 10 ( 𝜋 𝑟 2 0 ) + 2 . 5 log 10 5 ∼ 28 mag arcsec -2 .  \n<!-- image -->  \nFigure 17. Top : Stacked UV surface brightness profile of 200 galaxies with 𝑀 ∗ ∼ 10 9 M ⊙ at 𝑧 ≃ 6 in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. The gray squares and shaded region show the median surface brightness profile and 1𝜎 scatter. We fit the profile with the Sérsic profile and test the impact of fixing 𝑛 to unity (an exponential profile) and applying a surface brightness limit of 28 mag arcsec -2 . The dashed lines show the best-fit profiles and the colored circles indicate the corresponding 𝑅 UV eff derived. Bottom : A comparison of the fiducial 𝑅 UV eff with the one measured using a non-parametric area-based approach and the ones with different levels of AGN (point source) contamination. The area-based approach results in slightly smaller galaxy sizes due to the clumping of stellar distributions. An ambitious assumption of 50 percent contamination from a point source can result in about a factor of two smaller galaxy sizes. But they are still larger than the observational constraints.  \n<!-- image -->  \nprofile with the Sérsic model. We find that the free fit is already biased by the central bright component and underpredicts the surface brightness at the wing, yielding a small Sérsic index 𝑛 ∼ 0 . 6 and 𝑅 eff smaller than the true value. Fixing 𝑛 = 1 results in larger fitted 𝑅 eff that is surprisingly more consistent with the true value, but the fitting overpredicts the surface brightness at the wing. If we limit the fitting to where the flux is above 28 mag arcsec -2 , the best-fit 𝑅 eff is further biased to smaller values and is about a factor of two smaller than the true 𝑅 eff . Similar trends are found for galaxies at other stellar mass scales as well. This highlights the potential uncertainties in measuring galaxy sizes by fitting the surface brightness distribution. Optimistically speaking, the simulation predictions and observations  \ncould be reconciled by a combination of selection effects and fitting biases associated with the surface brightness limit.  \nClumpiness : In addition, the clumpy nature of high-redshift galaxies can also bias the size measurements to the separation of individual star-forming complexes or young star clusters (e.g. Vanzella et al. 2017, 2023; Chen et al. 2023; Claeyssens et al. 2023). Individual clump sizes will be substantially smaller than the sizes reported by classical dark matter halo finders. As discussed in e.g. Curtis-Lake et al. (2016); Ma et al. (2018a), an alternative way to measure size is by calculating the effective area 𝑆 eff (not necessarily continuous) enclosing half of the total surface brightness. The effective radius is then 𝑅 eff ≡ √︁ 𝑆 eff / 𝜋 . This non-parametric approach will be more sensitive to the sizes of individual stellar clumps and less biased by their separations. We compare this with the fiducial method in the bottom panel of Figure 17. We find that the size measured with this approach can be slightly smaller than our fiducial approach, but not enough to explain the discrepancy between simulations and observations. However, we cannot rule out the possibility that /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc underpredicts the clumpiness of the stellar distribution. The inability to resolve the multiphase ISM structure of galaxies in the IllustrisTNG model can result in different sub-kpc scale morphology of galaxies. Future studies using zoom-in simulations employing explicit ISM models are required to better understand the impact of clumpiness on galaxy size measurements.  \nUnresolved sources : AGN acts as a point source in galaxy photometry. In the observational sample, Morishita et al. (2024) used the [OIII]-H 𝛽 ratios and the mass-extinction diagram to demonstrate that their compact sources are not AGN or have little light contribution from AGN. From simulations, we calculate the bolometric luminosities of the SMBH in simulations as  \n𝐿 bol = ( 1 -𝜖 fb ) 𝜖 r / 𝑀 BH 𝑐 2 , (9)  \nwhere 𝜖 fb = 0 . 1is the feedback efficiency in the high-accretion quasar mode (Weinberger et al. 2017) and 𝜖 r is the radiative efficiency, assumed to be 0 . 2 in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc (Kannan et al. 2022a). We then convert this to the rest-frame UV luminosity using the bolometric corrections in Shen et al. (2020a). We neglect dust extinction here and the AGN luminosity should be interpreted as an upper limit.  \nAs previously shown in section 4.2, SMBH accretion rates are found higher in massive compact galaxies. However, even neglecting dust extinction, we find that the UV luminosities of AGN in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 are at least one order of magnitude smaller than the UV luminosities of galaxies at all stellar masses. If the AGN luminosity is subdominant compared to the galaxy luminosity, it should not have any impact on the half-light radius of the galaxy. But this statement strongly depends on the physics model for SMBH growth and feedback employed in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc, and is subjected to great uncertainties. The actual impact of unresolved AGN contamination could be significant. For example, in the case study of GN-z11, the two-component (a point source plus an extended component) fitting leads to about three times larger effective radius of the galaxy (Tacchella et al. 2023b).  \nWe can estimate the potential impact of unresolved sources by computing the radius containing 45 (25) percent of the total light, which corresponds to the contamination of a point source at the centre of the galaxy with 10 (50) percent of the total galaxy light. We show the comparison in the bottom panel of Figure 17. We find that 50 percent of contamination can result in about a factor of two smaller galaxy sizes across the whole stellar mass range. However, even with this ambitious assumption of point source contamination, this is not enough to explain the discrepancy with observations.', '7 CONCLUSIONS': 'In this paper, we provide a thorough analysis of galaxy intrinsic and apparent sizes at 𝑧 ≳ 6 based on the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations. We study the correlation of galaxy size with various galaxy properties and the large-scale environments they reside in. We investigate the physical mechanisms that determine the angular momentum transport in galaxies and drive the compaction in massive galaxies. Through comparison with observations, we find a discrepancy between the galaxy size predictions of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc and observations at the low-mass/faint end. Our key findings can be summarized below.  \n- · Intrinsic size-mass relation: The intrinsic half-mass radius increases (decreases) with galaxy stellar mass in low-mass (massive) galaxies. The breaking point of the two trends is around 10 8 M ⊙ and is insensitive to redshift. The redshift evolution at the low-mass end is mild. The most massive galaxies in the simulation volume continuously become more compact from 𝑧 ≃ 12 to 𝑧 ≃ 6 and eventually reach 𝑟 ∗ 1 / 2 ≲ 200 -300 pc. In the ( 95 . 5 cMpc ) 3 volume of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1, four of these compact massive galaxies at 𝑧 ≃ 6 are found to be quenched. Although the limitation of numerical resolution of the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 simulation could make the simulation miss a subpopulation of extremely compact galaxies, the median size mass relations derived from our simulations show good agreement with the high-resolution /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc simulation using identical physics and other cosmological simulations with different physics input and resolution.\n- · Correlation with galaxy properties: More compact galaxies tend to have lower halo spin at fixed stellar mass, which qualitatively agrees with the classical disk formation theory. In addition, more compact galaxies have lower sSFR, stronger AGN activity, lower gas abundance, and larger gas disk fraction. This suggests stronger feedback regulation after a period of centralized star-formation and rapid gas depletion, likely driven by the gravitational instability in gas disks. By tracking the evolution of the main progenitors of galaxies selected at 𝑧 ≃ 6, we find that compact (quenched) galaxies typically undergo three phases of size evolution, the initial growth of mass and size driven by mergers and accretion, the rapid gas depletion and compaction phase, and quenching after reaching a characteristic stellar surface density.\n- · Correlation with environmental properties: Galaxies in denser environments tend to have more extended stellar distributions. The evolution track of massive galaxies all show an early phase of size growth with frequent minor or major mergers, and later more compact galaxies have quieter merger histories.\n- · Comparison to analytical models: A deeper analysis of the size-spin correlation shows that the sizes of massive galaxies agree roughly with the prediction of the analytical disk formation theory. Though the most massive compact galaxy population experiences substantial loss in specific angular momentum. However, the sizes of low-mass galaxies barely correlate with the halo spin and are more consistent with a spherical shell model with galaxy sizes set by the shock radius of cold gas inflow and feedback-driven outflow. Therefore, the hump-like feature in the intrinsic size-mass relation can be explained by a merger-driven, spherical shock model for size growth of low-mass galaxies and a disk instability-driven compaction track of massive galaxies.\n- · Gas angular momentum profile: We investigate the specific angular momentum profiles of dark matter, gas, and stars of simulated galaxies, We find a clear decoupling of gas and dark matter at  \n- 0 . 1 𝑅 vir ≲ 𝑟 ≲ 0 . 3 𝑅 vir and the slow alignment of gas angular momentum with the central stellar disk. This phenomenon likely corresponds to the inflowing cold gas streams that roughly conserve angular momentum before mixing with the lower-angular momentum gas accreted earlier at radii below the impact parameter. The stellar specific angular momentum is about 3 times smaller than that of gas, which is a consequence of both feedback-driven outflow and disk instability. Low-mass galaxies have poorer spin alignment with the overall spin of the galaxy and this is consistent with their sizes better determined by the spherical shell model.\n- · Discrepancy with observations: We compute the apparent sizes of galaxies based on mock images generated using the Monte-Carlo radiative transfer method. The negative slope of the intrinsic size-mass relation is overturned by dust attenuation and we find a mild positive correlation between the UV effective radius and galaxy stellar mass as found in observations. A combination of projection effects, tracer effects, dust attenuation, and band selection can lead to order unity differences between the measured galaxy size and the intrinsic one, with non-trivial stellar mass dependence. The median apparent sizes of galaxies in rest-frame UV are at least three times larger than the observed ones at 10 7 M ⊙ ≲ 𝑀 ∗ ≲ 10 9 M ⊙ or -20 ≲ 𝑀 UV ≲ -15 at 6 ≲ 𝑧 ≲ 10. Such a discrepancy is robust against numerical uncertainties and exists in many other cosmological simulations with various physics inputs, simulation techniques, and numeric resolutions. Other than physical reasons, the selection effects and fitting biases associated with the surface brightness limit of existing galaxy surveys could be the major observational factors driving this discrepancy.  \nThe size-mass and size-luminosity relations during the epoch of reionization in the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations show distinct features compared to those observed at lower redshifts. These relations reveal a two-phase behavior where the sizes of low-mass galaxies are significantly influenced by the strength of feedback-driven outflows, whereas the sizes of more massive galaxies are affected by internal disk instability. These findings are especially relevant in the era of JWST, providing an opportunity to enhance our understanding of galaxy formation through observations of galaxy sizes or morphology in general at high redshifts.', 'ACKNOWLEDGEMENTS': 'The computations in this paper were run on the Engaging cluster at Massachusetts Institute of Technology (MIT). MV acknowledges support through the National Aeronautics and Space Administration (NASA) Astrophysics Theory Program (ATP) 19-ATP19-0019, 19-ATP19-0020, 19-ATP19-0167, and the National Science Foundation (NSF) grants AST-1814053, AST-1814259, AST-1909831, AST-2007355, and AST-2107724. EG acknowledges support from the CANON Foundation Europe through the Canon Fellowship program during part of the work presented in this paper.  \nSoftware citations:  \n- · N/u.pc/m.pc/p.pc/y.pc: Harris et al. (2020)\n- · S/c.pc/i.pc/p.pc/y.pc: Virtanen et al. (2020)\n- · A/s.pc/t.pc/r.pc/o.pc/p.pc/y.pc: Astropy Collaboration et al. (2013, 2018, 2022)\n- · M/a.pc/t.pc/p.pc/l.pc/o.pc/t.pc/l.pc/i.pc/b.pc: Hunter (2007)\n- · S/w.pc/i.pc/f.pc/t.pc/s.pc/i.pc/m.pc/i.pc/o.pc: Borrow & Borrisov (2020); Borrow & Kelly (2021)\n- · A/r.pc/e.pc/p.pc/o.pc-/r.pc/t.pc: Springel (2010); Kannan et al. (2019); Weinberger et al. (2020)', 'DATA AVAILABILITY': 'The data underlying this paper can be shared upon request to the corresponding author of this paper. The /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulation data is publicly available at https://www.thesan-project.com/. The postprocessed data and analysis scripts are stored on the Engaging cluster at MIT.', 'REFERENCES': "Vogelsberger M., Genel S., Sĳacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031  \nVogelsberger M., et al., 2014a, MNRAS, 444, 1518  \nVogelsberger M., et al., 2014b, Nature, 509, 177  \nVogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020a, Nature Re- views Physics, 2, 42  \nVogelsberger M., et al., 2020b, MNRAS, 492, 5167  \nWechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52  \nWeinberger R., et al., 2017, MNRAS, 465, 3291  \nWeinberger R., Springel V., Pakmor R., 2020, ApJS, 248, 32  \nWetzel A. R., Nagai D., 2015, ApJ, 808, 40  \nWhite S. D. M., 1984, ApJ, 286, 38  \nWiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99  \nWindhorst R. A., et al., 2023, AJ, 165, 13  \nWise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, MNRAS, 442, 2560  \nWolf J., Martinez G. D., Bullock J. S., Kaplinghat M., Geha M., Muñoz R. R., Simon J. D., Avedo F. F., 2010, MNRAS, 406, 1220  \nWu X., Davé R., Tacchella S., Lotz J., 2020, MNRAS, 494, 5636  \nXu H., Wise J. H., Norman M. L., Ahn K., O'Shea B. W., 2016, ApJ, 833, 84 Yajima H., Nagamine K., Zhu Q., Khochfar S., Dalla Vecchia C., 2017, ApJ, 846, 30  \nYajima H., et al., 2022, MNRAS, 509, 4037 Yang L., et al., 2022, ApJ, 938, L17 Yeh J. Y. C., et al., 2023, MNRAS, 520, 2757 Zolotov A., et al., 2015, MNRAS, 450, 2327  \nvan den Bosch F. C., Abel T., Croft R. A. C., Hernquist L., White S. D. M., 2002, ApJ, 576, 21  \nvan der Wel A., et al., 2014, ApJ, 788, 28", 'A1 Gravitational softening': "The choice of gravitational softening imposes an effective limit on the accelerations that may be adequately reproduced in the system. The density and circular velocity profiles can deviate from the true solution when the specific radial acceleration of particles, 𝑉 2 c ( 𝑟 )/ 𝑟 , drops below a characteristic value of 𝑎 𝜖 ∼ 0 . 5 𝐺𝑉 2 vir / 𝜖 (Power et al. 2003). Therefore, an arbitrarily large 𝜖 can break the convergence of mass density distribution at larger radii. This typically happens at a radius of an order unity factor times 𝜖 .  \nHowever, for typical total mass density profiles of the halo (e.g. NFW, Hernquist), a well-defined maximum value exists for radial acceleration. If 𝜖 is chosen such that 𝑎 𝜖 is always larger than this maximum value, the softening will not further affect the convergence of the mass density profile. For example, for the NFW profile, this condition is satisfied when  \n𝜖 ≲ ln ( 1 + 𝑐 ) -𝑐 /( 1 + 𝑐 ) 𝑐 2 𝑅 vir ≃ 0 . 05 𝑅 vir , (A1)  \nwhere we assume 𝑐 = 4 for typical high-redshift haloes. We obtain the halo mass limit where this effect becomes important as  \n𝑀 𝜖 vir = 2 . 4 × 10 7 M ⊙ × GLYPH<18> 𝜖 0 . 3 pkpc GLYPH<19> 3 . (A2)  \nFor the 𝜖 of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 and /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc, all the haloes relevant to our analysis in the main text are safely above this threshold.  \nThe softening itself acts as a way of 'smoothing' mass distribution and 'puffing up' structures. A useful reference scale is 𝜖 sp ≃ 2 . 8 𝜖 , which is the radius of the compact support given the cubic spline kernel chosen by the simulations. This is where the gravitational  \nacceleration becomes purely Newtonian. However, decreasing softening length does not lead to decreasing (and converging) galaxy sizes given the same particle mass resolution. First, an arbitrarily small softening length will inflate the integration error given the same time-stepping choice (Power et al. 2003). Secondly, small softening lengths exacerbate the segregation of stars and dark matter particles in halo centers and can counterintuitively increase galaxy sizes as softening length reduced (Ludlow et al. 2019a, 2020, 2023), which will be discussed in sections below. In simulations with sufficient mass resolution (Ludlow et al. 2020, 2023), galaxy sizes can become smaller than the canonical 𝜖 sp limit and converge in relatively massive galaxies ( 𝑀 ∗ ≳ 100 𝑚 b , where 𝑚 b is the baryonic mass resolution) when 𝜖 is chosen to be slightly smaller than 𝑟 conv (introduced below).", 'A2 Two-body relaxation': "When a finite number of particles is used to represent a system, individual particle accelerations will inevitably deviate from the mean-field value when particles pass close to each other. Even when orbits are integrated with perfect accuracy, these 'collisions' lead to changes of the order of unity in energy on the relaxation time scale (Binney & Tremaine 1987)  \n𝑡 relax ( 𝑟 ) ≃ 𝑡 circ ( 𝑟 ) 𝑁 ( < 𝑟 ) ln ( 𝑟 / 𝜖 ) , (A3)  \nwhere 𝑡 circ ≡ 2 𝜋𝑟 / 𝑉 c ( 𝑟 ) is the circular orbital time scale. This is more sensitive to the particle number than the gravitational softening length. Power et al. (2003) suggests the following convergence criterion for halo density profile against two-body relaxation  \n𝑡 relax ( 𝑟 ) 𝑡 circ ( 𝑅 vir ) ≡ √ 200 8 𝑁 ( < 𝑟 ) ln 𝑁 ( < 𝑟 ) GLYPH<18> ¯ 𝜌 𝜌 crit GLYPH<19> -1 / 2 ≳ 𝜅 , (A4)  \ngiven optimized choices of softening, force accuracy, and time stepping. 𝜅 is an order unity constant that has been studied in detail in Ludlow et al. (2019b), and a simpler expression of 𝑟 conv was found as 𝑟 conv ≃ 0 . 055 𝑙 /( 1 + 𝑧 ) with little halo mass dependence, where 𝑙 is the comoving mean interparticle separation in cosmological simulations. Taking the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 simulation as an example, 𝑟 conv takes the value of about 0 . 36 pkpc at 𝑧 ≃ 6.", 'A3 Energy transfer between particle species': 'In addition to two-body relaxation, close encounters between collisionless particles in simulations can also transfer energy between different species with different particle mass (Ludlow et al. 2019a, 2020, 2023). In the /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc simulations, the dark matter particle mass is set five times larger than the initial mass of gas cells (and roughly speaking the typical mass of stellar particles). Energy equipartition of the system will tend to transfer energy from dark matter to baryonic particles and spuriously increase the size of galaxies. This can lead to a flattening in galaxy size mass relation at a scale similar to 𝑟 conv (see Figure 7 in Ludlow et al. 2020) independent of 𝜖 .  \nLudlow et al. (2023) provided a thorough study of the typical time scale for spurious heating. By comparing this time scale with the age of the Universe, they derived convergence criteria for a range of galaxy properties including 𝑟 ∗ 1 / 2 . We choose 𝑁 lim 200 = 30 as defined in Table 2 in Ludlow et al. (2023) as a conservative estimate. We also scale it according to the change in Hubble time 1 / 𝐻 ( 𝑧 ) , which is a proxy for the age of the Universe at 𝑧 . We obtain the minimum halo mass scale that simulated galaxies are safe against spurious  \nFigure A1. The relation between galaxy sSFR and stellar mass in simulations with different resolutions. The solid line shows the median relation while the dashed lines show 1𝜎 scatter. The vertical lines show the mass limits corresponding to 50 times the baryonic mass resolution of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. The sSFR converges above the mass limit we choose. At lower masses, sSFR can be artificially and systematically suppressed in low-resolution simulations.  \n<!-- image -->  \nFigure A2. Size-mass relation in simulations with different resolutions at 𝑧 ≃ 6. The plotting style is the same as Figure A1. The two blue arrows indicate the convergence radius and the spurious heating mass limit of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1, respectively. We find good convergence of the size-mass relation at the resolution level of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1.  \n<!-- image -->  \nheating and convert it to stellar mass with a 𝑀 ∗ / 𝑀 vir ratio of 0 . 03 (conservatively large), which will be referred to as 𝑀 lim spur . This mass scale is marked with arrows in Figure A2 below and Figure 2 in the main text.', 'A4 Convergence of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 and /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc simulations': 'In Figure A1, we show the galaxy sSFR versus stellar mass relation at 𝑧 ≃ 6 in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 and two /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc simulations (L8-N512 and L4-N512). The numerical details of these simulations are listed in Table 1. We find that the sSFR converges at around 50 times the  \nFigure B1. Size-mass relation of central versus satellite galaxies at 𝑧 ≃ 6 in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 (colored distribution) and /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc (solid points). Although the satellites are slightly smaller than the centrals, the difference is too small to cause any notable shift in the median size of central galaxies versus all galaxies.  \n<!-- image -->  \nbaryonic mass resolution of the simulation. Below this mass, the sSFRs of galaxies are systematically suppressed in low-resolution simulations and a population of artificially quenched galaxies can show up. Therefore, in the galaxy classification step described in Section 3.2, we have specifically removed this population from our analysis. In Figure A2, we show the galaxy size-mass relation at 𝑧 ≃ 6 in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1 and /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-/h.pc/r.pc simulations. We find that both the median and 1𝜎 scatter of the relation converge already at the resolution level of /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc-1. The sizes and stellar masses of simulated galaxies are safely above the spurious heating regime defined by 𝑟 conv and 𝑀 lim spur calculated above.', 'APPENDIX B: CENTRAL VERSUS SATELLITE GALAXIES': 'In Figure B1, we show the size-mass relation of central versus satellite galaxies at 𝑧 ≃ 6 in /t.pc/h.pc/e.pc/s.pc/a.pc/n.pc. At fixed stellar mass, the satellites are slightly smaller than the centrals. However, this difference is small enough and it results in almost no differences in the median size of central galaxies and that of all galaxies. We also note that the small region around 𝑀 ∗ ∼ 10 7 . 5 -10 8 . 5 M ⊙ , 𝑟 ∗ 1 / 2 ∼ 0 . 5 pkpc is purely dominated by satellite galaxies. This explains the low-mass compact galaxy population at the same position in Figure 6, which resides in dense environments and is likely affected by the tidal fields and ram pressure of their hosts.', 'APPENDIX C: SIZE-MASS RELATION IN DIFFERENT ENVIRONMENTS': 'In Section 4.3 in the main text, we present the galaxy size-mass relation color-coded by the total matter overdensity and redshift of reionization, smoothed over 250 ckpc scale. Here in Figure C1, we show the results when the smoothing length is increased to the protocluster scale (1 cMpc). Regarding overdensity, we find the same trend as in Figure 6. Galaxies in underdense regions tend to be more  \n<!-- image -->  \nFigure C1. Galaxy intrinsic size-mass relation at 𝑧 ≃ 5 . 5 color-coded by the matter overdensity (top), and redshift of reionization (bottom) in the 1 cMpc environment of galaxies. This is the same as Figure 6 except for larger smoothing length when computing overdensity of 𝑧 reion .  \n<!-- image -->  \ncompact at fixed stellar mass. However, the dependence on 𝑧 reion is reversed and is completely controlled by the overdensity. The redshift of reionization on cMpc scales is strongly correlated with the amount of radiation sources in the environment and thus mass overdensity. Denser regions are reionized earlier. As summarized in the main text, the difference in the trends identified at 250 ckpc versus 1 cMpc scales reflects the transition of large-scale reionization (which depends more on overdensity and the ionizing photon budget in the environment) to small-scale reionization (where local ionizing sources and gas abundance start to be more relevant).  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}