Recent JWST observations with superb angular resolution have revealed the existence of clumpy galaxies at high redshift through the detection of restframe optical emission lines. We use the FirstLight simulation to study the properties of subgalactic clumps that are bright in the O III 5007 line with flux greater than 10SUP18SUP erg sSUP1SUP cmSUP2SUP to be detected by JWST. For 62 simulated galaxies that have stellar masses of 0.56  10SUP10SUP M SUBSUB at z  5 we find clumps in 1828 snapshots in the redshift range z  9.55.5. The clumps are identified by the surface density of the star formation rate SFR. About onetenth of the snapshots show the existence of clumpy systems with two or more components. Most of the clumps are formed by mergers and can be characterized by their ages central clumps dominated by stellar populations older than 50 Myr and offcentered clumps dominated by younger stellar populations with specific SFRs of 50 GyrSUP1SUP. The latter type of young clumps is formed from gas debris in the tidal tails of major mergers with baryonic mass ratios of 1  q lt 4. The mergerinduced clumps are shortlived and merge within a dynamical time of several tens of million years. The number density of the clumpy systems is estimated to be 10SUP5SUP cMpcSUP3SUP which is large enough to be detected in recent JWST surveys.
['High-redshift galaxies', 'Hydrodynamical simulations', 'James Webb Space Telescope', 'Galaxy formation', 'Reionization', '734', '767', '2291', '595', '1383', 'Astrophysics - Astrophysics of Galaxies']
A Mergerdriven Scenario for Clumpy Galaxy Formation in the Epoch of Reionization Physical Properties of Clumps in the FirstLight Simulation
{'A merger-driven scenario for clumpy galaxy formation in the epoch of reionization: Physical properties of clumps in the FirstLight simulation': 'Yurina Nakazato, 1 Daniel Ceverino, 2, 3 and Naoki Yoshida 1, 4, 5  \n1 Department of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan  \n- 2 Universidad Autonoma de Madrid, Ciudad Universitaria de Cantoblanco, E-28049 Madrid, Spain  \n3 CIAFF, Facultad de Ciencias, Universidad Autonoma de Madrid, E-28049 Madrid, Spain  \n- 4 Kavli Institute for the Physics and Mathematics of the Universe (WPI), UT Institute for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan\n- 5 Research Center for the Early Universe, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan', 'ABSTRACT': 'Recent JWST observations with superb angular resolution have revealed the existence of clumpy galaxies at high redshift through the detection of rest-frame optical emission lines. We use the FirstLight simulation to study the properties of (sub-)galactic clumps that are bright in [ Oiii ] 5007 ˚ A line with flux greater than ∼ 10 -18 erg s -1 cm -2 , to be detected by JWST. For 62 simulated galaxies that have stellar masses of (0 . 5 -6) × 10 10 M ⊙ at z = 5, we find clumps in 1828 snapshots in the redshift range z = 9 . 5 -5 . 5. The clumps are identified by the surface density of star formation rate. About one-tenth of the snapshots show the existence of clumpy systems with two or more components. Most of the clumps are formed by mergers and can be characterized by their ages; central clumps dominated by stellar populations older than 50 Myr, and off-centered clumps dominated by younger stellar populations with specific star formation rates of ∼ 50 Gyr -1 . The latter type of young clumps is formed from gas debris in the tidal tails of major mergers with baryonic mass ratios of 1 ≤ q < 4. The merger-induced clumps are short-lived, and merge within a dynamical time of several tens million years. The number density of the clumpy systems is estimated to be ∼ 10 -5 cMpc -3 , which is large enough to be detected in recent JWST surveys.', '1. INTRODUCTION': "Recent observations by the James Webb Space Telescope (JWST) have significantly advanced our understanding of galaxy formation and evolution in the early Universe. A number of highz galaxies have been observed so far, and the record of the most distant galaxies has been renewed (e.g., Curtis-Lake et al. 2023; Wang et al. 2023; Arrabal Haro et al. 2023; Harikane et al. 2023a; Carniani et al. 2024). For galaxies at z ≳ 9 . 5, strong rest-frame optical emission lines such as [ Oiii ] 5007 ˚ A are outside the range of JWST NIRSpec. Additionally, because of the vast distance, these galaxies appear typically to be compact single-component systems (e.g., Ono et al. 2023; Tacchella et al. 2023). Thus, it is difficult to obtain detailed information on kinematics and internal structure.  \nCorresponding author: Yurina Nakazato  \nyurina.nakazato@phys.s.u-tokyo.ac.jp  \nGalaxies at z = 6 -9 are promising targets for detailed studies on clump formation and evolution using emission lines. NIRCam has a resolution of 0.03'/pixel for long wavelength (2.4-5.0 µ m), corresponding to ∼ 150 pc/pixel at z = 8 1 . NIRSpec Integral Field Unit (IFU) has a pixel scale of 0.05' 2 , corresponding to ∼ 250 pc/pixel at z = 8. Gravitational lensing with magnification µ increases the effective spatial resolution by √ µ , allowing detailed studies of internal structure of galaxies (e.g., Hsiao et al. 2023a; Bradaˇc et al. 2023; ' Alvarez-M'arquez et al. 2023; Jones et al. 2023; Hashimoto et al. 2023; Morishita et al. 2023; Meˇstri'c et al. 2022; Claeyssens et al. 2023; Fujimoto et al. 2024).  \n- 2 NIRSpec IFU elements are 0.1'/pixel on the sky, but the combination of dithering and drizzle weighting allows sub-sampling of the detector pixels with 0.05' (Jones et al. 2023; Hashimoto et al. 2023).  \nEspecially, recent JWST photometry observations have reported a large number of galaxies at z > 6 that contain star-forming clumps and /or elongated structures with a scale of ∼ 100 pc in the rest-UV and optical wavelengths. (e.g., Hainline et al. 2023; Chen et al. 2023; Huertas-Company et al. 2024; Tacchella et al. 2023; Baker et al. 2023). A prominent example is JADESGN-189.18051+62.18047, which consists of five clumps, and four of them are located within a small region of ∼ 7 kpc.  \nSome NIRSpec IFU observations also discovered merging systems of star-forming galaxies and core regions of proto-clusters at z > 6 with the strong detection of [ Oiii ] 5007 ˚ A(e.g., Hashimoto et al. 2023; Jones et al. 2023; Arribas et al. 2023; Venturi et al. 2024). Hashimoto et al. (2023) and Jones et al. (2023) found over four components within a ∼ 11 kpc regions. While these observations suggest a rapid assembly of high redshift galaxies, it remains largely unknown how the compact clumpy systems emerge in the early epochs.  \nClump formation at lower redshift z = 0 -4 have been widely investigated both by observations (Genzel et al. 2011; Tadaki et al. 2014; Guo et al. 2015; Elmegreen et al. 2013) and theoretical studies (Oklopˇci'c et al. 2017; Ceverino et al. 2010; Bournaud & Elmegreen 2009; Agertz et al. 2009; Romeo et al. 2010; Romeo & Agertz 2014). They have proposed two major mechanisms for clump formation: violent disk instability (VDI) and mergers.  \nThe former is driven by, for example, cold gas accretion flowing along the cosmic web (Dekel et al. 2009; Ceverino et al. 2010; Dekel & Krumholz 2013; Rodriguez-Gomez et al. 2016), and the local disk instability is regulated by the Toomre parameter Q (Toomre 1964). This is originally defined as Q ≡ σ r κ/πG Σ for an infinitesimally thin gas disk, where σ r , κ, G, Σ are the radial velocity dispersion, the epicyclic frequency, the gravitational constant, and gas surface density. Clumps formed in massive disk galaxies have been studied by cosmological simulations (Agertz et al. 2009; Ceverino et al. 2010; Mandelker et al. 2014, 2017; Buck et al. 2017; Ceverino et al. 2023; Inoue et al. 2016; Genel et al. 2012) as well as isolated disk models (Tamburello et al. 2017; Faure et al. 2021; Bournaud et al. 2014; Bournaud & Elmegreen 2009; Hopkins et al. 2012). Observations also identified clumps induced by VDI with significant rotations (e.g., Genzel et al. 2006, 2008). The fate of formed clumps (survival and disruption) has often been discussed from theory and observations (Dekel et al. 2009; Renaud et al. 2014; Elmegreen et al. 2007).  \nThe later competitive mechanism is mergers. Galaxy mergers cause significant gas density perturbations and  \neven promote star formation (e.g., Di Matteo et al. 2008), and thus may be another promising formation path for clumps. Such clumpy galaxies via mergers are often discussed with a major-merger rate or a pair fraction. Numerical simulations suggest that the galaxy merger rate may increase at higher redshifts (e.g., Genel et al. 2008, 2009; Fakhouri & Ma 2009; Hopkins et al. 2010; Rodriguez-Gomez et al. 2015).  \nSince observations cannot measure the rate directly, they instead quantify pair fractions (e.g., Tasca et al. 2014; Ribeiro et al. 2017; L'opez-Sanjuan et al. 2013), which has been statistically investigated at lower redshift z ≲ 6 with a large samples by HST. Ribeiro et al. (2017) argued that the obtained pair fraction infers a major merger fraction of ∼ 20 % over the redshift range 2 < z < 6.  \nIn addition to a pair fraction, a clumpy fraction f clumpy has also been investigated to evaluate the contribution of merger and VDI (e.g., Guo et al. 2015; Shibuya et al. 2016; Ribeiro et al. 2017). Ribeiro et al. (2017) found that a clumpy fraction increases from 35 to 50 percent from z = 2-6, while Shibuya et al. (2016) reported the opposite trend and suggested that VDI is the main origin of clumps. The clumpy fraction may also depend on stellar mass (Guo et al. 2015), making it difficult to understand the redshift evolution.  \nWhile clumpy galaxies at z > 7 are newly observed by JWSTin multiwavelength and their stellar mass and size are statistically investigated (e.g., Meˇstri'c et al. 2022; Claeyssens et al. 2023), the redshift evolution of clumpy fractions and constraints on the formation scenario are not well understood. For numerical simulations, largescale cosmological simulations such as Illustris TNG (Pillepich et al. 2018), THESAN (Shen et al. 2024), FLARES (Roper et al. 2022), BLUETIDES (Marshall et al. 2022), and SIMBA (Wu et al. 2020) have a minimum resolution of a few hundred pc, making it hard to resolve internal structures of galaxies at high redshift.  \nZoom-in simulations can achieve better resolution, and clumpy structures are also reported (Barrow et al. 2017; Ma et al. 2018a; Katz et al. 2019; Arata et al. 2020; Gelli et al. 2020, 2021; Kohandel et al. 2020; Kimm et al. 2016; Calura et al. 2022). However, their samples are often limited to only a few galaxies, and there are few theoretical studies at z > 6 that qualitatively and statistically discuss the formation of clumps and their fates, their physical properties, and the clumpy fractions.  \nIn this paper, we use the FirstLight simulation suite, which achieves a very high resolution of 17 pc and provides a large sample size to allow robust statistical analysis. We use outputs for 62 simulated galaxies with 4092 snapshots in z = 6 -9. Our main objective is to study  \nthe properties of galaxy clumps at the epoch of reionization, and their formation and evolution.  \nIn our previous study (Nakazato et al. 2023), we devised a physical model of H ii regions and calculated the luminosity of emission lines for the same FirstLight galaxy samples. The resulting relation between the star formation rate and [ Oiii ] luminosity is consistent with recent observations (e.g., Hashimoto et al. 2018; Tamura et al. 2019; Harikane et al. 2020). We use the same model to characterize clumpy galaxies that are bright in rest-frame optical emission lines such as [ Oiii ] 5007 ˚ A.  \nThe rest of the present paper is organized as follows. In Section 2, we introduce the numerical simulation and describe the method to calculate line luminosities and define clumps. In Section 3, we discuss in detail the formation and evolution of a few clumpy systems. Summary and conclusions are given in Section 4. Throughout this paper, we assume a set of parameters in ΛCDM cosmology with Ω m = 0 . 307 , Ω b = 0 . 048 , h = 0 . 678 and σ 8 = 0 . 823 from Planck13 results (Planck Collaboration et al. 2014).", '2.1. Zoom-in Simulation': 'We use the cosmological zoom-in simulation suite FirstLight (Ceverino et al. 2017). The simulations are performed using the art code (Kravtsov et al. 1997; Kravtsov 2003; Ceverino & Klypin 2009), which follows gravitational N -body dynamics and Eulerian hydrodynamics with Adaptive Mesh Refinement (AMR). The parent simulation has a box size of 40 comoving h -1 Mpc on a side. We identified the 62 most massive halos in the parent dark matter-only simulation. The halos have maximum rotation velocities with V max > 178 kms -1 at z = 5. We then performed zoom-in simulations in these regions with an effective resolution of 4096 3 and a maximum resolution of 17-32 proper pc. The dark matter particle mass is m DM = 8 × 10 4 M ⊙ and the minimum stellar particle mass is 10 3 M ⊙ . We have compared the stellar-to-halo mass relation with other simulations with the same setup but a higher resolution of m DM = 10 4 M ⊙ and confirmed that the convergence is satisfactory (Ceverino et al. 2023, 2024). A number of baryonic processes in star formation and feedback are incorporated in a sub-grid manner. We describe them briefly as follows.  \nThe zoom-in simulations follow radiative cooling by atomic hydrogen and helium, metal ions and atoms, and molecular hydrogen, and photoionization heating by a time-dependent UV background with partial selfshielding (Haardt & Madau 1996).  \nStellar particles are created in gas cells according to the criteria for the hydrogen number density of n H > 1 cm -3 and temperature of T < 10 4 K. We use a stochastic model with the probability that scales with the gas free-fall time (Ceverino & Klypin 2009). Our simulations incorporate a constant heating rate for 40 Myr after star formation as thermal stellar feedback driven by stellar winds and supernovae (SNe). The code also implements radiative feedback from massive young stars as a non-thermal pressure term (Ceverino et al. 2014). The kinetic feedback model also includes the injection of momentum coming from the (unresolved) expansion of gaseous shells from supernovae and stellar winds (Ostriker & Shetty 2011). More details can be found in Ceverino et al. (2017), Ceverino et al. (2014), Ceverino et al. (2010), and Ceverino & Klypin (2009).  \nOur simulation implements a pressure floor to prevent artificial fragmentation (Truelove et al. 1997; Machacek et al. 2001; Robertson & Kravtsov 2008; Agertz et al. 2009; Rosdahl et al. 2015). In order to resolve Jeans length at least seven cells ( N J = 7), we set the artificial pressure floor at  \nP floor = Gρ 2 gas N 2 J ∆ x 2 min πγ , (1)  \nwhere ∆ x min and γ = 5 / 3 are the minimum cell size and the specific heat ratio for monoatomic gas. The pressure in the Euler equation is replaced by an effective pressure that is set to P floor when it would have a lower value. The threshold number of cells (i.e., N J = 7) is beyond the standard criteria ( N J = 4, Truelove et al. 1997) and this shows good convergence in cosmological simulations of clumpy galaxies (Ceverino et al. 2010). Note that this implementation prevents artificial fragmentation but it also prevents the growth of structures below ∼ 100 pc scales. Previous simulations with the same implementation and resolution have already investigated low-redshift clumps (Mandelker et al. 2014, 2017; Inoue et al. 2016; Inoue & Yoshida 2019; Ceverino et al. 2023). More details are provided in Appendix B.  \nOur 62 galaxy samples have stellar masses greater than ∼ 10 10 M ⊙ at z = 5. The maximum resolution in the zoom-in hydro simulations, particularly valid in dense clumps, enables the simulations to resolve gas densities of ∼ 10 3 cm -3 with temperatures of ∼ 300 K. We store a total of 66 snapshots between z = 9 . 5 to z = 5 . 5 for each galaxy, with a spacing of the cosmic expansion parameter ∆ a = 0 . 001, corresponding to 7 -10 Myr 3 . The frequent output timing is sufficient to follow dy-  \ning galaxy mergers that occur on the tidal timescale of 10-100 Myrs (Renaud et al. 2013).', '2.2. Line luminosity Calculation': "We calculate the luminosities of emission lines from H ii regions for the simulated galaxies and clumps. We adopt essentially the same model as Nakazato et al. (2023) and assume that each stellar particle is surrounded by a spherical H ii region. We generate a table of a set of emission lines using cloudy (Ferland et al. 2013). As an input stellar SEDs, we adopt a single star model of BPASS (Eldridge et al. 2017). According to Eldridge et al. (2017), these SEDs are similar to other ones such as Starburst 99 (Leitherer et al. 1999), GALAXEV (Bruzual & Charlot 2003), and Maraston (2005) at the age of < 1 Gyrs. Note that BPASS binary star model can enhance the number of ionizing photons by a factor of up to ∼ 10 at the age of < 10 Myr, leading to the [O iii ] flux increase up to ∼ 1.2 times higher than that of the single model (Ceverino et al. 2021). The input parameters for cloudy are gas metallicity Z gas , ionization parameter U , and gas density in H ii regions n hii . The adopted values of the parameters are listed in Table 1. In our table, we normalize the emission line luminosity for each stellar particle ( L line ,i ) by H β luminosity as  \nL line ,i = (1 -f esc ) C line ( Z gas , U, n hii ) L caseB H β . (2)  \nHere f esc is the escape fraction and we assume f esc = 0 . 1, a value suggested by both radiation hydrodynamics simulations (Yajima et al. 2011; Kimm & Cen 2014; Wise et al. 2014; Trebitsch et al. 2017; Xu et al. 2016), and observations (e.g., Nakajima et al. 2020) for highz galaxies. Since individual H ii regions are not resolved in our simulations, we assume that the gas density of a H ii region is given by  \nn hii = { 100 cm -3 ( n gas , grid < 100 cm -3 ) n gas , grid ( n gas , grid ≥ 100 cm -3 ) , . (3)  \nwhere n gas , grid is the gas number density of a grid with a size of ∆ = 50 pc. Our motivation for using a variable H ii region density comes from the need to account for grid-to-grid density variations in our simulations. Setting a fixed H ii region density at a typical value of 100 cm -3 (Osterbrock & Ferland 2006; Hirschmann et al. 2017, 2022) could lead to unrealistic conditions where the grid density exceeds this value. To address this, we adopt a model with a variable H ii region density  \nif n gas , grid > 100 cm -34 . Our model results in a higher global-averaged H ii region density in highz ( z > 6) simulated galaxies compared to typical local H ii regions. This trend is also supported by recent JWST observations (e.g., Isobe et al. 2023; Fujimoto et al. 2022; Abdurro'uf et al. 2024).  \nWe look for the line ratio table C line with the three parameters Z gas , U, n hii which are closest to those of the target stellar particles in our simulation. We obtain the line luminosity for each grid as L line , grid = ∑ N i L line ,i , where N is the number of stellar particles that the grid contains.  \nTable 1. The parameters used to calculate the line luminosities with cloudy .  \nIn our previous study (Nakazato et al. 2023; Mushtaq et al. 2023), we find that the magnitude of dust attenuation at rest-frame 5007 ˚ A is approximately A 5007 = 0 . 08 , 0 . 24 , 0 . 7 for individual galaxies with M ∗ /M ⊙ = 10 8 , 10 9 , 10 10 , respectively. Thus, we expect that dust attenuation affects the line luminosity of a small clump by a factor of a few. We further discuss this issue of dust attenuation in Section 4.", '2.3. Clump Finder': "In order to compare the properties of the simulated clumpy galaxies with observations, we identify clumps in two-dimensional 'images' using projection along a random line of sight ( z -axis) of a cube with a fixed side length of 10 kpc. The size is larger than the main central galaxy, and is close to a field of view of NIRSpec IFU at z = 8. We assign masses of gas, stars, and dark matter to a uniform grid with ∆ = 50 pc by using a cloud-in-cell interpolation also employed in Mandelker et al. (2014, 2017) and Ceverino et al. (2023). The value is about three times the maximum AMR resolution. We first mark grids that contain young ( < 10 Myr) stars with their surface mass density (Σ M ∗ (young) ) greater than 10 8 . 5 M ⊙ kpc -2 , corresponding to the surface SFR density of 10 1 . 5 M ⊙ yr -1 kpc -2 . We then identify candidate groups by running a Friends-of-Friends finder with linking length ∆ to the marked grids. Groups with at least  \nFigure 1. The [ Oiii ] 5007 ˚ A luminosity vs. SFR for our 62 simulated galaxies at z = 9 (purple), 8 (blue), 7 (green), and 6 (yellow).  \n<!-- image -->  \n16 grids ( N grid ) are identified as clumps. The threshold grouping number of grids corresponds to the minimum radius of clump R min = √ N grid ∆ 2 /π = 113pc. The threshold surface stellar mass density and the minimum number of grids are chosen to ensure that the identified clumps are resolved by observations by NIRSpec IFU and NIRCam with the line flux of F [ Oiii ] ≳ 10 -18 erg s -1 cm -2 at z = 8. The details of threshold determination are described in the next section 2.3.1. We have also checked the effect of line of sight by changing inclination, and it does not change the number of clumps.", '2.3.1. Parameter determination': "We set three parameters for our clump finder; grid size (∆), threshold surface mass density of young stellar population (Σ M ∗ (young) ), and the minimum number of grids ( N grid ). We assume that a H ii region exists around each stellar particle. The grid size is set to be larger than the typical Stromgren radius for each stellar particle in our line emission model. The Stromgren radius R s is given by  \nR s = ( 3 Q 0 4 πn 2 HII α B ) 1 / 3 (4)  \n= 5 . 4 pc ( n HII 10 2 cm -3 ) -2 / 3 ( Q 0 5 × 10 49 s -1 ) 1 / 3 , (5)  \nwhere n HII is the gas density in H ii regions and Q 0 is the number of ionizing photons per unit time. The electron temperature is T e = 10 4 K, and a filling factor of ϵ = 1 is assumed. The value of Q 0 depends on stellar age, mass, and spectrum, which means the value varies to simulation resolutions and model SEDs. We  \nhave checked the distribution of Q 0 values for all stellar particles in our simulations and obtained the median (maximum) value of Q 0 = 5 × 10 50 (5 × 10 52 ) s -1 , and the corresponding Stromgren radius is R s ∼ 15 (54) pc for n HII = 100cm -3 . We thus set ∆ = 50 pc.  \nTo determine the threshold of surface young stellar mass density (Σ M ∗ (young) ), we consider the relationship between L [ Oiii ] 5007 ˚ A and star formation rate (SFR). This relationship is implied from pre-JWST observations that have shown the tight correlation between [O iii ] 88 µ m and SFR at z > 6 (e.g., Harikane et al. 2020; Moriwaki et al. 2018; Katz et al. 2019; Pallottini et al. 2022). There are still few statistical investigations of [O iii ] 5007 ˚ A emitters at z > 6 (Matthee et al. 2023; Meyer et al. 2024), but lower-redshift ( z ≲ 2) galaxies are already reported to have tight correlation between [O iii ] 5007 ˚ A and SFR by both observations and simulations (e.g., Villa-V'elez et al. 2021; Hirschmann et al. 2023). Figure 1 shows the [ Oiii ] 5007 ˚ Aluminosity against SFR for our galaxy samples. It is approximated as 5  \nL [ Oiii ] 5007 ˚ A ∼ 6 × 10 7 ( SFR M ⊙ yr -1 ) L ⊙ (6)  \n∼ 6 × ( M ∗ (young) M ⊙ ) L ⊙ , (7)  \nwhere M ∗ (young) is the mass of very young ( < 10 Myrs) stars, and the instantaneous SFR is estimated as ( M ∗ (young) / 10 7 ) M ⊙ yr -1 . The flux limit to detect a clump by [ Oiii ] 5007 ˚ A line is ∼ 3 × 10 -18 erg s -1 cm -2 , which corresponds to signal-to-noise ratio of S/N > 5 with an exposure time of ∼ 10 4 sec (Hashimoto et al. 2018; Jones et al. 2023; Arribas et al. 2023; Matthee et al. 2023). The corresponding luminosity is L [ Oiii ] 5007 ˚ A ∼ 6 × 10 7 L ⊙ at z = 8. Substituting this into eq.(7), we obtain M ∗ , (young) ∼ 10 7 M ⊙ for a 'detectable' clump.  \nWe set the minimum radius to be R clump ∼ 100 pc, which can be resolved by recent JWST NIRCam/NIRSpec IFU observations (e.g. Chen et al. 2023; Hainline et al. 2023; Hashimoto et al. 2023; Harikane et al. 2024). We therefore set the threshold grid grouping number as N grid = 16, corresponding to R clump = √ N grid ∆ 2 /π = 113pc, where one grid area is ∆ 2 =  \n5 Our linear fitted line for each redshift is  \nlog 10 L [ Oiii ] [ L ⊙ ] = 1 . 16 × log 10 SFR + 7 . 89 ( z = 9)  \nlog 10 L [ Oiii ] [ L ⊙ ] = 1 . 16 × log 10 SFR + 7 . 74 ( z = 8)  \nlog 10 L [ Oiii ] [ L ⊙ ] = 0 . 87 × log 10 SFR + 7 . 97 ( z = 7)  \nlog 10 L [ Oiii ] [ L ⊙ ] = 0 . 97 × log 10 SFR + 7 . 67 ( z = 6) .  \n50 2 pc 2 . The minimum clump size is around seven times larger than the minimum AMR resolution of 17 pc. We emphasize that our identified clumps correspond to unlensed observed clumps. Strong lensing can enhance the resolution up to the ∼ 10 pc level (e.g. Adamo et al. 2024; Bradley et al. 2024; Mowla et al. 2024).  \nFor this minimum clumps to have 'detectable' luminosity, i.e., to have young stellar mass of M ∗ (young) ≳ 10 7 M ⊙ , each grid should have a surface stellar mass density of Σ M ∗ (young) ≥ 10 7 / 16 M ⊙ grid -1 = 10 8 . 5 M ⊙ kpc -2 . This value also corresponds to the surface SFR density of Σ(SFR) ≥ 1 / 16 M ⊙ yr -1 grid -1 = 10 1 . 5 M ⊙ yr -1 kpc -2 .", '3.1. Projection of physical properties for simulated clumpy galaxies': 'Figure 2 shows three examples of simulated galaxies. Two, four, and three clumps are identified in these outputs as bright components. The clumps trace the distribution of young stars but also that of dense gas with n gas ≳ 100 cm -3 . Tail features can also be seen in [ Oiii ] emission, similar to the recent galaxies observed by JWST such as MACS1149-JD1 ( z = 9 . 1, Bradaˇc et al. 2023) and B14-65666 ( z = 7 . 15, Sugahara et al. 2024). The distribution of the stars older than 100 Myrs appears more diffuse than the young stellar population. The underlying dark matter distribution shows signatures of mergers, which will be explained in detail in Section 3.3. Some clumps (e.g., clump C in Figure 2) are located off-center. The clumps are separated with ∼ 0 . 1 -4 kpc distances, similarly to the clumps at z ≃ 6 -8 studied by Chen et al. (2023), with typical separations of ∼ 0.3 - 4.3 kpc. Such close clumps are also found in other JWST observations at z > 6; the JADES (Hainline et al. 2023), EIGER (Matthee et al. 2023) and CANUCS (Asada et al. 2024) surveys. The stellar age distributions in the right-most panels show a clear difference in the distribution of young and old components. Interestingly, similar features are found by JWST MIRI observations (e.g., Colina et al. 2023).  \nThe panels in the middle row focus on the galaxy FL957, which contains the highest number of clumps at z > 6. To see if those clumps are gravitationally bound or not, we have estimated the virial parameter (Bertoldi & McKee 1992), α vir ≡ 5 σ 2 clump R c /GM , where σ clump and R c are velocity dispersion within a clump and clump radius. All clumps have α ∼ 0 . 2, which indicates that these clumps are gravitationally bound.  \nThe four-clump feature motivates us to compare it with the recent observation by Hashimoto et al. (2023), who spectroscopically studied the core of the most dis-  \ntant protocluster at z = 7 . 88, A2744-z7p9OD, with JWST NIRSpec integral field spectroscopy. The core region contains four clumps which are detected by [ Oiii ] 5007 ˚ A in a region of ∼ 11 × 11 kpc. Table 2 summarizes the clump properties. We find that the stellar mass, SFR, and [ Oiii ] luminosities are similar to the observed clumps (see Table 1 in Hashimoto et al. (2023)). These properties are typical for clumps younger than 50 Myrs identified in the multiple systems, as shown in Table 3, which is discussed in the next section. We discuss that such properties are different from clumps older than 50 Myr in Section 3.2. In Sections 3.3 and 3.4, we show the formation and evolution of multi-clump systems.', '3.2. Distributions of clump properties': 'Our clump finder identifies 2254 clumps in total. Of these, 1503 clumps are single systems, i.e., galaxies with one clump, often appearing as proto-bulge components. The remaining 751 clumps are in 325 clumpy systems. Table 3 summarizes the physical properties of the clumps. Figure 3 shows the distribution of the clump age weighted by stellar mass for single and multiclump objects. The distribution for single systems has a peak at ∼ 200 Myr, while the clumps in clumpy systems have a peak at ∼ 90 Myr. Based on the age distribution, we can separate them into two types of clumps: clumps younger than 50 Myr and clumps older than 50 Myr, which are indicated by blue and red histograms, respectively. Table 3 also shows the properties of clumps younger than 50 Myr and older than 50 Myr, respectively. We see that old clumps have similar stellar mass, age, and clump size to the clumps identified in single systems. Young clumps are smaller and their stellar mass is ten times lower than older clumps. Figure 4 shows the distributions of three properties for different systems: single systems with one clump ( N clump = 1, black solid line), and young and old clumps in clumpy systems ( N clump ≧ 2, red and blue). The probability distributions are normalized to the total number of clumps in each of the three systems.  \nThe left panel of Figure 4 refers to the gas fraction ( f gas ≡ M gas /M baryon ) of individual clumps. Single systems have a distribution in the range of f gas ∼ 0 . 4 -0 . 95 and have the peak at roughly f gas ∼ 0 . 85. While clumps younger than 50 Myrs concentrate on the value of 0.95, clumps older than 50 Myrs have almost the same distribution as single systems.  \nThe middle panel of Figure 4 refers to the surface mass ratio of baryon and dark matter. Single systems have a peak at Σ baryon / Σ DM ∼ 10, and the clumps older than 50 Myr have almost the same distribution and peak value. The shape of the distribution of clumps younger  \nFigure 2. Three simulated clumpy galaxies at z = 9 . 00 (top, ID=FL977), z = 7 . 70 (middle, ID=FL957), and z = 7 . 20 (bottom, ID=FL958). Each column from left to right represents the averaged gas density, surface stellar mass density of young ( < 10 Myr) population, [ Oiii ] 5007 ˚ Aline flux, surface mass density of old ( > 100 Myr) stellar population and dark matter, and mass-weighted stellar age, respectively. Each panel has a side length of 10 proper kpc and the same size of projected depth.  \n<!-- image -->  \nTable 2. Summary of properties of FL957 at z = 7 . 70, which contains four clumps within 10 × 10 kpc region. Note that M ∗ , (young) refers to stellar masses younger than 10 Myrs and M ∗ , (all) are stellar masses for all ages. A clump radius R c in the 5th column is obtained as R c = √ ∆ 2 N grid /π .  \nTable 3. Summary of physical properties of galaxies we identify in single systems and clumpy systems. The clumps in clumpy systems are further categorized into two groups: those with ages younger than 50 Myr and older than 50 Myr. The last row of the table lists the average values for both systems. We see that old clumps ( > 50 Myr) in clumpy systems have similar properties to ones in single systems, as discussed in section 3.2.  \nFigure 3. The distribution of mass-weighted clump ages. Black and dashed-purple histograms represent single systems (detected one central clump, i.e., proto-bulge) and clumpy systems, respectively. The histograms are normalized by the total number of clumps for each system. We divide the clumpy system into clumps younger than 50 Myrs (blue) and those older than 50 Myrs (red).  \n<!-- image -->  \nthan 50 Myr is similar to the old ones, but the peak is shifted to a 0.5 dex larger value. This implies that over 50% of young clumps are formed off-center of their host dark matter halos, i.e., in a baryonic-rich environment.  \nThe right panel of Figure 4 refers to the sSFR of individual clumps. The peak of the distribution for single systems is at sSFR ∼ 6 Gyr -1 , while the peaks for old and young clumps are at ∼ 10 and 55 Gyr -1 , respectively. The sSFR values of single systems and old ( > 50Myr) clumps in clumpy systems are similar to that of main-sequence galaxies at z ∼ 6, sSFR ∼ 6 Gyr -1 (Popesso et al. 2023). However, young clumps in clumpy systems have over 9 times larger sSFR than that of mainsequence galaxies, and such large values are also seen in observed merging galaxies at z > 6 (e.g., Sugahara et al. 2024).  \nFrom the three physical properties in Figure 4, we see that clumps older than 50 Myr tend to have the same properties as single systems. On the other hand, clumps younger than 50 Myr are located in baryonic-rich environments. They are also gas-rich and in the process of bursty star formation. This indicates a different formation path for young clumps.', '3.3. Clump formation induced by merger': "In this section, we study the formation of clumps during a merger. The top panels of Figure 5 show the time evolution of dark matter distribution around FL957, and the bottom panels represent gas density distribution in the same central region of 10 kpc. At z = 8 . 01, we see two dark matter halos, one in the center of the panel and the other one in the right bottom. The two galaxies experience the first pericenter passage at z = 7 . 85, which enhances the tidal compression of the surrounding gas and forms dense clumps. After the first passage ( z = 7 . 70), we identify four clumps and this number is larger than that of merging halos. These clumps are formed off-center from the halos. Tail structures coming from clumps are formed by tidal effects during the galaxy interaction 6 . The formation of clumps triggered by galaxy merger has been investigated by isolated simulations which initially set an Antennae-like galaxy merger (e.g., Renaud et al. 2013, 2014, 2015). These previous studies have found that the first two pericenter passages increase compressive (curl-free) turbulence and lead to clumpy morphologies and bursty star formation with a timescale of 10-30 Myr. Their results are consistent with our findings from cosmological simulations at high redshift. Unlike the isolated simulations, which set initial parameters for galaxies and their orbits to reproduce the morphology and kinematics of merging galaxies (e.g., Renaud et al. 2013, 2015; Fensch et al. 2017; Maji et al. 2017), our cosmological simulations naturally demonstrate the formation of clumps due to mergers occurring at high redshift. Our simulations can also investigate statistics by connecting to a major merger rate density (see section 3.5).  \nWe estimate the mass of a perturber that generates clumps based on the 0th-order approach of a spherical collapse model for a homogeneous sphere (Steidel et al. 1998; Overzier 2016). A merger-induced clump is formed if the mass is  \nM clump = ¯ ρδ × V clump ≳ M min , (8)  \nwhere M min is the minimum baryon mass of the observable clump we define. Here, ¯ ρ, δ, V clump are the mean gas density around the progenitor, typical overdensity induced by the merger, and volume of the formed clump. Assuming that the perturber compresses the surrounding gas contained in a volume with its size ∼ R per , the formed clump mass is re-written as  \nM clump = 2 × 10 8 M ⊙ ( ¯ ρδ 250 cm -3 )( R per 200 pc ) 3 . (9)  \n<!-- image -->  \n<!-- image -->  \nFigure 4. The probability distribution of the gas fraction, baryon-to-dark matter mass ratio, and specific star formation rate for our clump samples. The histograms with black solid lines are for single systems, and we color the histograms for clumps depending on their age. The distributions are normalized to the total number of clumps in each of the colored types.  \n<!-- image -->  \nFigure 5. Time evolution of dark matter density (top) and gas density distribution (bottom) from z = 8 . 01 to z = 7 . 33. The side and projected lengths are 10 kpc for both projections. The black circles represent the location of the identified clumps. We set t = 0 at z = 7 . 70, when the galaxy consists of four luminous clumps.  \n<!-- image -->  \nWe normalize the relevant quantities by typical values found in our simulations (Table 3). For example, ¯ ρδ ∼ 250 cm -3 . The minimum mass of a clump that emits [ Oiii ] lines bright enough to be observed is derived based on the young stellar mass criteria M ∗ (young) , min = 10 7 M ⊙ in section 2.3.1,  \nM min ∼ M ∗ (young) , min 1 -f gas = 10 8 M ⊙ ( 1 -f gas 0 . 1 ) -1 , (10)  \nwhere f gas is the gas fraction defined in section 3.2. Young clumps have f gas ∼ 0 . 9 (Figure 4). Substituting equations (9) and (10) into eq.(8), we obtain the characteristic size of the baryonic 'wake' induced by  \nmerger encounter (perturber) as R per ≳ 150 -200 pc. The mass ratio q ≡ M center /M per largely determines the subsequent morphological evolution and the mode of triggered star formation. We assume the central one is more massive, i.e., M center ≥ M per . Considering that the central object in the system has a typical size of R center ≃ 200 -230 pc as listed in Table 3 and assuming that both objects have similar gas densities, we estimate the mass ratio to be q = ( R center /R per ) 3 = 1 -3 . 6. We thus argue that bright clumpy galaxies in optical emission lines at z > 6 can be formed by a major merger.  \n3.4. Clump fate  \nFigure 6. Time evolution of clump number as a function of redshift for the three galaxies shown in Figure 2. The inverted triangles represent the redshifts in Figure 2. N clump = 0 represents that clumps are too faint to be identified by the clump finder.  \n<!-- image -->  \nFigure 5 also follows the subsequent evolution of the clumpy system. A further 18 Myr after the formation of the four clumps, Clump B and C merge with clump A and D resulting in two clumps, which finally merge into a single system 26 Myrs later.  \nFigure 6 shows the redshift evolution of the number of clumps for the three galaxies shown in Figure 2. Clumpy structures are observed temporarily and finally merge into the center of the system 30-40 Myr after the multiple systems are formed, whose evolution is seen in other galaxy samples. The typical timescale is the crossing time during the merger:  \nt cross = d v -v center ∼ 50 Myr ( d 3 kpc )( v -v center 60 km / s ) -1 , (11)  \nwhere d and v center are the separation and the velocity of the central object. The derived crossing timescale is consistent with the merger timescale in Figure 5, and roughly the same order of the orbital (dynamical) timescale of this system in a 10 kpc region, t orb ∼ t dyn ∼ 25 Myr. Therefore, we see that clumpy systems which are bright in rest-frame optical emission lines are temporary structures and finally merge into one massive galaxy within ∼ 50 Myr.  \nNotice that the two galaxy systems (FL957, FL977) in Figure 6 are identified as N clump = 0 at z ∼ 9. This implies that these galaxies are too faint to be identified as an [ Oiii ] bright object at z ∼ 9. As noted in Section 2, our clump identification is based on surface SFR density and the threshold value is 30 M ⊙ yr -1 kpc -2 (cor-  \nding to 1 M ⊙ yr -1 per clump). However, a lensing effect enables us to identify clumpy features even at z ≳ 9 as the recent observations report triply lensed galaxies at z = 10 . 17 (Hsiao et al. 2023a,b).  \nSince our simulations end at z = 5, we cannot follow the further evolution to the present day. Some galaxies form clumps via disk instabilities at z < 6, outside the redshift range investigated here. Such clumps still exist at the final snapshot, z = 5, but most probably they migrate to the center or they are disrupted due to feedback. Due to our limited resolution, we cannot exclude the possibility that the densest regions of the clumps formed proto-globular clusters that may survive all these events. Future modeling will be required to follow their evolution.", '3.5. number abundance of clumpy galaxies': "Figure 7 represents the fraction of clump number at each redshift. Around 90% of galaxies at z = 8 -6 are observed as a single system ( N clump = 1). From Table 2 and 3, we see that both observable clumpy galaxies tend to have stellar masses of M ∗ ≳ 10 9 M ⊙ in their systems, even though our clump definition sets a minimum clump mass of M ∗ ≥ M ∗ (young) , min = 10 7 M ⊙ . Our sample galaxies at z ≥ 8 . 5 have small stellar masses of M ∗ ≲ 10 9 M ⊙ , resulting in that all of our samples at 8 . 5 ≤ z ≤ 9 are identified as a single system. At these highredshifts, the corresponding UV magnitude is M UV ≳ -19 . 7 7 . Galaxies at z ≲ 8 are bright ( M UV < - 19.7) and multiple clumps are identified, which are formed by mergers. This trend is interestingly consistent with the recent JWST observations of Chen et al. (2023), who find that luminous z ≃ 6 -8 galaxies with M UV < -20 . 7 tend to have 3-4 clumps and fainter galaxies with M UV > -20 . 7 tend to be a single system. The fraction of galaxies with three or more clumps ( N clump ≧ 3) is ∼ 0.5 - 1%.  \nWe do not find any clumpy systems formed by violent disk instability (VDI) through cold accretion until z ∼ 6. This is because VDI requires a massive disk with a large disk mass fraction (Dekel et al. 2009; Ceverino et al. 2010). We identify some disk-like galaxies at z = 8 -9, but the stellar masses are ≲ 5 × 10 9 M ⊙ and have smoothed disk. Specifically, two galaxies with  \n7 We use the scaling relation between stellar mass and UV magnitude obtained from Ceverino et al. (2017, 2019) as follows,  \nM UV = log 10 M ∗ -6 + M ∗ ∗ (for z = 6 -7) , (12)  \nα ∗ M UV = log 10 M ∗ -6 α ∗ + M ∗ ∗ -1 (for z = 8 -9) , (13)  \nwhere α ∗ = -0 . 394 ± 0 . 002 and M ∗ ∗ = -12 . 13 ± 0 . 03.  \nFigure 7. The fraction of each clump number at z ∼ 9 (purple), z ∼ 8 (blue), z ∼ 7 (green), and z ∼ 6 (yellow). The clump number N clump is referred to as the number of clumps in the galaxy system within 10 kpc × 10 kpc region.  \n<!-- image -->  \nfive clumps are seen in a massive system; the total stellar mass in the field is M ∗ ∼ 10 10 M ⊙ at z ∼ 5 . 5. The clumps in such a massive system are expected to be formed via violent disk instability as same as lowerredshift galaxies (Mandelker et al. 2014, 2017; Ceverino et al. 2023; Inoue et al. 2016).  \nFigure 8 shows the number density of clumpy galaxies. As described in Section 2.1, our study employs zoom-in simulations and we need to correct for incompleteness in the number count of clumpy galaxies. We use the results of a larger boxsize N -body simulation (Klypin et al. 2011) (refer to Appendix A for details). The number density of single systems is (1 -3) × 10 -4 cMpc -3 at z = 6 -9. We identify all the simulated galaxies at z = 9 as single systems and the number density is consistent with the UVLF at M UV ∼ -20 (Harikane et al. 2023b). At z = 6 -8, we estimate the number density to be (1 -5) × 10 -5 cMpc -3 , and this value is enough to be observed by recent JWST surveys. The EIGER survey (Kashino et al. 2023) conducted deep JWST/NIRCam wide-field slitless spectroscopic observations and identified 117 [ Oiii ] emitters at z = 5 . 33 -6 . 93 (Matthee et al. 2023). They found that ∼ 1% of their observed [ Oiii ] emitters are closely separated within < 1' (corresponding to ≲ 6 kpc at z ∼ 6). They also derived the [ Oiii ] luminosity function at z ∼ 6. We integrate it within the luminosity range for their samples and obtain the number density of their observed [ Oiii ] emitters, n [ Oiii ] = 10 -2 . 8 cMpc -3 . Thus, we can estimate the number density of clumpy systems of their observed [ Oiii ] emitters as ∼ 0 . 01 × n [ Oiii ] = 1 . 6 × 10 -5 cMpc -3 , consistent with our theoretical abundance of clumpy galaxies at z = 6 -8.  \nFigure 8. The number density of clumpy galaxies as a function of redshift. The gray and purple plots represent single systems and clumpy systems, respectively. The error bars are calculated from Poisson errors, √ N/V box , where N is the number of clumpy galaxies and V box is the comoving volume of the simulation box.  \n<!-- image -->  \nredshift  \nWe can also derive the merger rate density from the clumpy systems within a 10 kpc × 10 kpc region following Mo et al. (2010). Under the assumption that these pairs will merge on a time scale τ mrg , the merger rate density is related to the clumpy factor f clumpy ≡ (# of clumpy systems with N clump ≧ 2)/(# of total systems), according to  \n˙ n mrg = 1 ⟨ N clump ⟩ f clumpy n gal τ mrg = 6 . 3 × 10 -5 cMpc -3 Gyr -1 × ( f clumpy 0 . 1 )( n gal 10 -4 cMpc -3 )( τ mrg 80 Myr ) -1 , (14)  \nwhere ⟨ N clump ⟩ is the averaged clump number among clumpy systems. The factor 1/ ⟨ N clump ⟩ takes into account the fact that a multi-clump system ends in a single merger. Figure 7 shows that ∼ 80 - 90% of clumpy systems consist of two clumps, and we thus adopt ⟨ N clump ⟩ = 2. We obtain the clumpy factor from Figure 8 as f clumpy ∼ 0 . 1 -0 . 2 and see that the value does not change very much at z = 6 -9. For the comoving number density of galaxies population n gal , we obtain n gal ∼ 10 -4 cMpc -3 with M ∗ ≳ 10 9 M ⊙ 8 from Figure 8. The value also does not change largely within a factor of 1-2. We adopt the crossing timescale (eq. (11)) as the merger timescale τ mrg . In our case, we take  \ninto account a galaxy pair located within a 10 kpc × 10 kpc region, i.e., the separation of d ∼ 5 kpc, which gives τ mrg ∼ 80 Myr. The merging timescale is also seen in the top panels of Figure 5. Adopting this merger as a typical case, we obtain the merger rate density of ˙ n mrg = 6 . 3 × 10 -5 cMpc -3 Gyr -1 . This value is roughly consistent with the results from Millennium Simulation at z > 6 (fig.2 in Kitzbichler & White 2008). We also compare our merger rate density with that of Illustris simulations. Rodriguez-Gomez et al. (2015) derive the major merger rate ∼ 1 Gyr -1 for the descendant stellar masses of ∼ 10 9 M ⊙ at z = 6 -9. From the stellar mass function at z > 6 of Illustris simulations (Genel et al. 2014), the corresponding number density of the descendant mass is ∼ 10 -4 cMpc -3 . Therefore, the major merger rate density is ∼ 10 -4 cMpc -3 Gyr -1 , which is in good agreement with our values. These consistencies infer that clumpy galaxies, which are bright in rest-frame optical emission lines, can be linked to major mergers.", '4. DISCUSSION & SUMMARY': "We have studied formation of clumpy galaxies in the epoch of reionization observable in rest-frame optical emission lines. The clumps we identify are bright in [ Oiii ] 5007 ˚ Aand have radii of ∼ 150 - 200 pc, which can be resolved by JWST. We have shown that the observable large clumps are formed by major mergers. The clumps are categorized into two types: proto-bulges dominated by stellar populations older than 50 Myr and off-centered clumps dominated by younger stellar populations. The latter type of clump is formed from gas debris in tidal tails induced by the mergers. The merger-driven clumps tend to have a high gas fraction of over 90%, trigger bursty star formation with sSFR ∼ 20 -100 Gyr -1 , and produce young stellar components. The clumpy galaxies account for ∼ 10% of the population, and the fraction remains roughly constant from z = 9 to z = 6. We also find that the clumpy systems are short-lived, ending in a merger with a companion(s) within a few tens of Myrs. The number density of clumpy systems is estimated to be (1 -5) × 10 -5 cMpc -3 , which is large enough to be observed in recent JWST surveys (Matthee et al. 2023).  \nFormation of clumpy galaxies at z ≲ 4 has been studied in detail (Mandelker et al. 2014, 2017; Buck et al. 2017; Ceverino et al. 2023; Inoue et al. 2016). Especially, Mandelker et al. (2014) study clump formation in 29 galaxy samples at z = 1 -4. They define clumps as baryonic overdensities in three dimensions. They find that 70% of the clumps are formed in situ by violent disk instability while the remaining 30% are formed ex situ  \nby minor mergers. In the present paper, we study clump formation at z = 6 -9 using high-resolution simulations. We define the clumps as [ Oiii ] bright objects based on the threshold surface star formation rate so that they can be observed by JWST with a reasonable exposure time. We find that most of the observable clumps are formed by mergers of galaxies with roughly the same baryonic mass. A sufficiently large 'perturber' comparable to the central object with size ∼ 200 pc can generate gas wakes and overdensities to form large, bright clumps. Interestingly, violent disk instability (VDI) is less likely to occur in low-mass galaxies because the disk mass itself is small (Dekel et al. 2009). In the VDI case, even when it occurs, clumps contain only 1-2 % of the mass of the central object (Ceverino et al. 2010), rendering them too faint to be detected through current observations. Our result suggests a change of the mode of clump formation from the predominantly merger-driven scenario in the early universe to VDI at low redshifts at z ≲ 4 -5.  \nThere are still some caveats in our study. Firstly, we define clumps by adopting a constant threshold of surface star formation rate density. The constant value is motivated to identify observable clumps, but at the same time, we might miss clumps in low-mass galaxies at z ∼ 9. Such unidentified low-mass clumps might correspond to proto-globular clusters. Recent JWST photometric observations in gravitational lens fields have found young stellar cluster candidates with M ∗ = 10 5 -10 6 M ⊙ at z = 6 -10 (Vanzella et al. 2022, 2023; Adamo et al. 2024; Bradley et al. 2024; Mowla et al. 2024). The star clusters are also very small in size with effective radii of less than 1 pc (Adamo et al. 2024), which our simulations cannot resolve. Some cosmological simulations focus on the formation of dwarf galaxies and stellar clusters at z > 6 and have sub-pc scale resolution (Ma et al. 2018b,c; Calura et al. 2022; Garcia et al. 2023; Sugimura et al. 2024), allowing one to study the formation of smaller-scale structures.  \nSecond, our clump identification is based on restframe optical emission lines. Statistics of highz ( z > 6) clumpy galaxies have been conducted from rest-frame UV observation by HST (e.g., Shibuya et al. 2016; Bowler et al. 2017, 2022; Meˇstri'c et al. 2022). Some samples have also been observed by ALMA, and they have offsets between FIR emission regions and UV emission regions (Carniani et al. 2017). Recent JWST multiband photometry or IFU observations enable us to perform multiwavelength studies of clumpy galaxies, from the rest-UV to NIR. The number of clumps is often reported to have spatial offsets between UV, optical, and NIR wavelengths (Colina et al. 2023; Bik et al. 2023;  \nKalita et al. 2024; Rodighiero et al. 2024), while Treu et al. (2023) find no dramatic morphology changes across multi-bands. We note that multi-component observations at high redshift ( z > 6) are more challenging than at lower redshift (e.g. z ∼ 2 -3, Meˇstri'c et al. 2022; Zanella et al. 2024), where galaxy morphologies are most of the times less disturbed and present clearer features.  \nOur study focuses on young giant clumps that are bright in [O iii ] because rest-optical emission lines are less affected by dust attenuation than UV emission. These lines trace dense star-forming regions that are gravitationally bound and exist as physical components (Claeyssens et al. 2023; Matthee et al. 2023). Clumps identified in UV-bands are expected to be heavily attenuated by dust in dense clouds of n gas ≥ 100 cm -3 (Figure 2). To detect UV clumps in the simulated galaxies and compare them with the current multiwavelength observations, we need to conduct post-processed threedimensional radiative transfer calculation (e.g., SKIRT, Baes et al. 2011; Camps & Baes 2020). We leave it as future studies (Nakazato et al. 2024 in prep. )  \nWe have found that a clumpy system contains different stellar/clump populations: old ( > 50 Myr) clumps at the center of dark matter halos, and young ( < 50 Myr) off-centered clumps. The latter type of clumps originates from dense gas clouds that are formed in tidal tails during mergers. Observationally, it would be interesting to identify old clumps with ages of a few hundred Myrs. As shown in Figure 2, there are spatial offsets between the young ( ≲ 50 Myr) and old ( ∼ 300 Myr) components. Those old components are not detected as a clump based on our criteria. Recent JWST observations show that some clumps are detected in rest-frame optical emission lines and continuum, while others are detected only in continuum (e.g., Hashimoto et al. 2023). Detailed radiative transfer calculations would also enable us to understand the cause of the spatial offsets between young stellar populations and dust continuum emitting regions as seen in recent observations at z > 6 by ALMA (Tamura et al. 2023).  \nEmission line flux ratios provide crucial information on clump properties such as the ISM electron density,  \nmetallicity, and ionization parameters. For instance, Yang et al. (2023) use a zoom-in simulation to show that the ratio of F [ Oiii ] 5007 ˚ A /F [ Oii ] 3727 ˚ A can be used to trace stellar populations younger than several Myr. The ratio is taken by some NIRSpec IFU/NIRISS observations to evaluate the ionization state and metalicity gradients in each pixel (Arribas et al. 2023; Saxena et al. 2024; Tripodi et al. 2024; Rodr'ıguez Del Pino et al. 2024; Wang et al. 2022; Venturi et al. 2024). Line ratios are also measured for z ∼ 4 lensed galaxies to determine and produce resolved dust-to-gas ratio maps (Birkin et al. 2023). Applying the same method to clumpy galaxies at higher redshift will enable us to characterize young and old clumps and to test the merger-driven scenario we propose here.", 'ACKNOWLEDGEMENTS': "We would like to thank the anonymous referee for very useful comments. We are grateful to thank Takatoshi Shibuya, Takuya Hashimoto, Yuma Sugahara, Javier ' Alvarez-M'arquez, Luis Colina, Luca Constantin, Carmen Blanco Prieto for the fruitful discussion. This work made use of v2.3 of the Binary Population and Spectral Synthesis (BPASS) models as described in Byrne et al. (2022) and Stanway & Eldridge (2018). The authors thankfully acknowledge the computer resources at MareNostrum and the technical support provided by the Barcelona Supercomputing Center (RES-AECT-2020-30019). Numerical analyses were carried out on the analysis servers at Center for Computational Astrophysics, National Astronomical Observatory of Japan. YN acknowledges funding from JSPS KAKENHI Grant Number 23KJ0728 and a JSR fellowship. DC is a RamonCajal Researcher and is supported by the Ministerio de Ciencia, Innovaci'on y Universidades (MICIU/FEDER) under research grant PID2021-122603NB-C21. YN and NYacknowledge support from JSPS International Leading Research 23K20035.", 'A. ABUNDANCE MATCHING': "FirstLight simulation provides a velocity-selected sample of massive galaxies with log V circ [km / s] ≥ 2 . 3 at z = 5. The output of the zoom-in simulations contains the galaxies and their progenitors but not other smaller galaxies in a cosmological volume. In order to estimate the comoving number density of low mass galaxies, we use the circular velocity function of Bolshoi N -body simulation, which has a larger volume of (250 h -1 cMpc) 3 (Klypin et al. 2011, 2016; Rodr'ıguez-Puebla et al. 2016). Figure 9 compares the velocity functions at z = 9 -6 obtained from FirstLight  \nFigure 9. Velocity function of halos at z =9(purple), 8(blue), 7(green), and 6(yellow). Solid lines are obtained from FirstLight 62 samples, and dashed lines are from the Bolshoi simulation using the same cosmological parameter (Klypin et al. 2011).  \n<!-- image -->  \n(solid lines) with those from the Bolshoi simulation (dashed lines). The latter can be well approximated by  \nn ( > V ) = AV -3 exp ( -[ V V 0 ] α ) , (A1)  \nwhere  \nA = 1 . 52 × 10 4 σ -3 / 4 8 ( z ) ( h -1 Mpc / kms -1 ) -3 , (A2)  \nα = 1 + 2 . 15 σ 4 / 3 8 ( z ) , (A3)  \nV 0 = 3300 σ 2 8 ( z ) 1 + 2 . 5 σ 2 8 ( z ) kms -1 , (A4)  \nas derived in Klypin et al. (2011).  \nWe reconstruct the halo and galaxy abundances in the following manner. We configure 20 velocity bins logarithmically scaled. In each bin, we calculate the ratio of the number density of the Bolshoi and FirstLight simulations, R i = n ( > V circ ) , Bolshoi /n ( > V circ ) , FL . For each snapshot subscripted i , we compute the circular velocity ( V circ ,i ) of the galaxy and count the total number of galaxies within a velocity bin. We then use the ratio R i as a weighting factor to derive the 'true' number of single (clumpy) systems at a given redshift as  \nN single(clumpy) ( z bin ) = #snap( z bin ) ∑ i Θ R i , (A5)  \nwhere Θ is set to 1 when snapshot i is single (clumpy) system and otherwise 0. Here #snap( z bin ) represents the number of snapshots in each redshift bin. Dividing the number by the volume of the FirstLight simulation, i.e., (40 cMpc /h ) 3 , we obtain the number density of the single or clumpy systems at each redshift in Figure 8.", 'B. NUMERICAL RESOLUTION AND PRESSURE FLOOR': 'In our simulations, we impose a pressure floor on the condensing gas in galaxies to prevent artificial fragmentation. It is known that fragmentation can occur if the local Jeans length λ J is not well resolved by computational grids with a cell size ∆ x min . Truelove et al. (1997) show that artificial fragmentation can be prevented if the Jeans length is always resolved by at least four resolution elements (i.e., N J ≡ λ J / ∆ x min = 4). Ceverino et al. (2010) suggest a more stringent criteria of N J = 7 by running long-term cosmological simulations.  \nFigure 10. Temperature-density phase diagram for a clumpy system FL958 at z = 7 . 2. The color bar shows the gas mass within each bin (∆ log 10 n gas , ∆log 10 T gas ) = (0 . 08 , 0 . 06). Gray dashed lines show the temperature threshold (horizontal) and density threshold (vertical) for star formation. The black dotted line shows the threshold of eq. (B7). We plot all of the gas cells inside a cube with a side length of 10 kpc.  \n<!-- image -->  \nIn our simulations, we set an effective pressure floor at  \nP floor = Gρ 2 gas N 2 J ∆ x 2 min πγ , (B6)  \nwith N J = 7. We assume a ratio of specific heats of γ = 5 / 3 for monoatomic gas and ρ gas is the gas mass density. The pressure in the Euler equation is replaced by eq. (B6) when its pure thermal pressure contribution has a lower value. Artificial fragmentation is effectively prevented by implementing the pressure floor, which is meant to model crudely non-thermal contributions such as local turbulence in the ISM.  \nFigure 10 shows a phase diagram for a snapshot of a clumpy system at z = 7 . 2 (see Figure 2). Gray dashed lines show the temperature threshold (horizontal, T crit = 10 4 K) and density threshold (vertical, n crit = 1cm -3 ) for star formation. The black dashed line in the left panel shows the following relationship calculated from eq. (B6),  \nT J = m p k B P floor ρ gas = 35K ( n gas 1 cm -3 ) ( N J 7 ) 2 ( ∆ x min 17 pc ) 2 , (B7)  \nwhere m p , k B and n gas are proton mass, Boltzmann constant, and gas number density, respectively. Note that there are gas cells below the black dashed line, on which the pressure floor is imposed in the simulation. The pressure floor has already been implemented in Ceverino et al. (2010) and detailed analysis on clumps at lower redshifts is given in some previous studies (Mandelker et al. 2014, 2017; Inoue et al. 2016; Inoue & Yoshida 2019; Ceverino et al. 2023) that use simulations with essentially the same set up.', 'REFERENCES': "Abdurro'uf, Larson, R. L., Coe, D., et al. 2024, arXiv e-prints, arXiv:2404.16201,  \ndoi: 10.48550/arXiv.2404.16201  \nAdamo, A., Bradley, L. D., Vanzella, E., et al. 2024, arXiv e-prints, arXiv:2401.03224,  \ndoi: 10.48550/arXiv.2401.03224  \nAgertz, O., Teyssier, R., & Moore, B. 2009, MNRAS, 397, L64, doi: 10.1111/j.1745-3933.2009.00685.x  \n' Alvarez-M'arquez, J., Colina, L., Crespo G'omez, A., et al.  \n2023, arXiv e-prints, arXiv:2309.06319,  \ndoi: 10.48550/arXiv.2309.06319  \nArata, S., Yajima, H., Nagamine, K., Abe, M., & Khochfar,  \nS. 2020, MNRAS, 498, 5541,  \ndoi: 10.1093/mnras/staa2809"}