Context. Feedback from stars in the form of radiation stellar winds and supernovae is crucial to regulating the star formation activity of galaxies. Dwarf galaxies are especially susceptible to these processes making them an ideal test bed for studying the effects of stellar feedback in detail. Recent numerical models have aimed to resolve the interstellar medium ISM in dwarf galaxies with a very high resolution of several solar masses. However when it comes to modeling the radiative feedback from stars many models opt for simplified approaches instead of explicitly solving radiative transfer RT because of the computational complexity involved. Aims. We introduce the Realistic ISM modeling in Galaxy Evolution and Lifecycles RIGEL model a novel framework to selfconsistently model the effects of stellar feedback in the multiphase ISM of dwarf galaxies with explicit RT on a starbystar basis. Methods. The RIGEL model integrates detailed implementations of feedback from individual massive stars into the stateoftheart radiationhydrodynamics code AREPORT. It forms individual massive stars from the resolved multiphase ISM by sampling the initial mass function and tracks their evolution individually. The lifetimes photon production rates massloss rates and wind velocities of these stars are determined by their initial masses and metallicities based on a library that incorporates a variety of stellar models. The RT equations are solved explicitly in seven spectral bins accounting for the infrared to He II ionizing bands using a momentbase scheme with the M1 closure relation. The thermochemistry model tracks the nonequilibrium H He chemistry as well as the equilibrium abundance of C I C II O I O II and CO in the irradiated ISM to capture the thermodynamics of all ISM phases from cold molecular gas to hot ionized gas. Results. We evaluated the performance of the RIGEL model using 1 MSUBSUB resolution simulations of isolated dwarf galaxies. We found that the star formation rate SFR and interstellar radiation field ISRF show strong positive correlations with the metallicity of the galaxy. Photoionization and photoheating can reduce the SFR by an order of magnitude by removing the available cold dense gas fuel for star formation. The presence of ISRF also significantly changes the thermal structure of the ISM. Radiative feedback occurs immediately after the birth of massive stars and rapidly disperses the molecular clouds within 1 Myr. As a consequence radiative feedback reduces the age spread of star clusters to less than 2 Myr prohibits the formation of massive star clusters and shapes the cluster initial mass function to a steep powerlaw form with a slope of   2. The massloading factor measured at z  1 kpc of the fiducial galaxy has a median of SUBMSUB  50 while turning off radiative feedback reduces this factor by an order of magnitude. Conclusions. We demonstrate that RIGEL effectively captures the nonlinear coupling of early radiative feedback and supernova feedback in the multiphase ISM of dwarf galaxies. This novel framework enables the utilization of a comprehensive stellar feedback and ISM model in cosmological simulations of dwarf galaxies and various galactic environments spanning a wide dynamic range in both space and time.
['hydrodynamics', 'radiative transfer', 'methods: numerical', 'ISM: general', 'galaxies: dwarf', 'galaxies: evolution', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Instrumentation and Methods for Astrophysics']
RIGEL Simulating dwarf galaxies at solar mass resolution with radiative transfer and feedback from individual massive stars
{'RIGEL: Simulating dwarf galaxies at solar mass resolution with radiative transfer and feedback from individual massive stars': '<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nYunwei Deng 1 , Hui Li 1 , Boyuan Liu 2 , Rahul Kannan 3 , Aaron Smith 4 , and Greg L. Bryan 5  \n- 1. Department of Astronomy, Tsinghua University, Haidian DS 100084, Beijing, China\n- e-mail: hliastro@tsinghua.edu.cn\n- 2. Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK\n- 3. Department of Physics and Astronomy, York University, 4700 Keele Street, Toronto, ON M3J 1P3, Canada\n- 4. Department of Physics, The University of Texas at Dallas, Richardson, Texas 75080, USA\n- 5. Department of Astronomy, Columbia University, New York, NY 10027, USA  \nReceived XXX XX, XXXX; accepted YYY YY, YYYY', 'ABSTRACT': 'Context. Feedback from stars in the form of radiation, stellar winds, and supernovae is crucial to regulating the star formation activity of galaxies. Dwarf galaxies are especially susceptible to these processes, making them an ideal test bed for studying the e ff ects of stellar feedback in detail. Recent numerical models have aimed to resolve the interstellar medium (ISM) in dwarf galaxies with a very high resolution of several solar masses. However, when it comes to modeling the radiative feedback from stars, many models opt for simplified approaches instead of explicitly solving radiative transfer (RT) because of the computational complexity involved.  \nAims. We introduce the Realistic ISM modeling in Galaxy Evolution and Lifecycles (RIGEL) model, a novel framework to selfconsistently model the e ff ects of stellar feedback in the multiphase ISM of dwarf galaxies with explicit RT on a star-by-star basis. Methods. The RIGEL model integrates detailed implementations of feedback from individual massive stars into the state-of-the-art radiation-hydrodynamics code, arepo -rt . It forms individual massive stars from the resolved multiphase ISM by sampling the initial mass function and tracks their evolution individually. The lifetimes, photon production rates, mass-loss rates, and wind velocities of these stars are determined by their initial masses and metallicities based on a library that incorporates a variety of stellar models. The RT equations are solved explicitly in seven spectral bins accounting for the infrared to He ii ionizing bands, using a momentbase scheme with the M1 closure relation. The thermochemistry model tracks the nonequilibrium H, He chemistry as well as the equilibrium abundance of C i , C ii , O i , O ii , and CO in the irradiated ISM to capture the thermodynamics of all ISM phases, from cold molecular gas to hot ionized gas.  \nResults. We evaluated the performance of the RIGEL model using 1 M ⊙ resolution simulations of isolated dwarf galaxies. We found that the star formation rate (SFR) and interstellar radiation field (ISRF) show strong positive correlations with the metallicity of the galaxy. Photoionization and photoheating can reduce the SFR by an order of magnitude by removing the available cold, dense gas fuel for star formation. The presence of ISRF also significantly changes the thermal structure of the ISM. Radiative feedback occurs immediately after the birth of massive stars and rapidly disperses the molecular clouds within 1 Myr. As a consequence, radiative feedback reduces the age spread of star clusters to less than 2 Myr, prohibits the formation of massive star clusters, and shapes the cluster initial mass function to a steep power-law form with a slope of ∼ -2. The mass-loading factor (measured at z = 1 kpc) of the fiducial galaxy has a median of η M ∼ 50, while turning o ff radiative feedback reduces this factor by an order of magnitude.  \nConclusions. We demonstrate that RIGEL e ff ectively captures the nonlinear coupling of early radiative feedback and supernova feedback in the multiphase ISM of dwarf galaxies. This novel framework enables the utilization of a comprehensive stellar feedback and ISM model in cosmological simulations of dwarf galaxies and various galactic environments spanning a wide dynamic range in both space and time.  \nKey words. galaxies: dwarf - galaxies: evolution - ISM: general - methods: numerical - hydrodynamics - radiative transfer', '1. Introduction': "Star formation on galactic scales is remarkably ine ffi cient (Kennicutt 1998; Bigiel et al. 2008; Sun et al. 2023), with only a small fraction of the total baryonic mass converted into stars in the Universe (Yang et al. 2003; Conroy & Wechsler 2009; Moster et al. 2013; Wechsler & Tinker 2018). This ine ffi ciency is particularly pronounced in dwarf galaxies ( M ⋆ < 10 9 M ⊙ ), which exhibit the highest dark-to-stellar mass ratios (e.g. de Blok & McGaugh 1997; Mateo 1998; Hayashi et al. 2020; Eftekhari et al. 2022) The low star formation e ffi ciency (SFE) in these low-mass systems is commonly attributed to various physical processes, including radiation, stellar winds, and supernovae (SNe) from massive stars, collectively referred to as 'stellar feedback.' They  \ndeposit mass, momentum, and energy into the cold interstellar medium (ISM), transforming it into a multiphase structure before potentially expelling it into the circumgalactic medium (CGM) or even the intergalactic medium (IGM) (e.g. McKee & Ostriker 1977; Dekel & Silk 1986; Hopkins et al. 2011, 2012b).  \nStellar feedback reduces galaxy-wide SFE by strongly shaping the evolution of small-scale star-forming regions. Feedback from young massive OB stars rapidly disrupts their natal giant molecular clouds (GMCs) and prevents further gravitational collapse and local star formation (Grudi'c et al. 2018; Li et al. 2019; Kim et al. 2021b). In particular, pre-SN (early) feedback mechanisms, especially radiation, play a crucial role as they contribute to the feedback budget that is comparable to or even larger than  \nthat of SN explosions (Agertz et al. 2013; Kannan et al. 2014; Geen et al. 2015; Hopkins et al. 2020; Je ff reson et al. 2021, see Schinnerer & Leroy 2024 for a review from the observational perspective). Moreover, recent observations indicate that radiative and stellar wind feedback disperses GMCs on very short timescales, well before SN events occur (Kruijssen et al. 2019; Chevance et al. 2022, 2023), making GMCs short-lived objects where stars and clusters formed within a very short time span of a few million years (Hartmann et al. 2012; Hollyhead et al. 2015; Kim et al. 2021a). On a larger scale, stellar feedback drives turbulence, stabilizing the galactic disk against runaway vertical collapse by balancing momentum injection with turbulent dissipation (Faucher-Giguère et al. 2013; Ostriker & Kim 2022). Consequently, the interplay between these local and global e ff ects regulates the galaxy-wide star formation (Dobbs et al. 2011; Hopkins et al. 2011; Semenov et al. 2017; Hayward & Hopkins 2017) and establishes an extended gas depletion timescale of 1-3 Gyr in star-forming galaxies (Bigiel et al. 2008; Leroy et al. 2008; Sun et al. 2023).  \nOver the past few decades, cosmological simulations have been able to reproduce the ine ffi ciency of cosmic star formation histories and various observed galaxy statistics, including the galaxy luminosity function, clustering, color bimodality, and scaling relations between galaxy mass, size, or metallicity hyperplanes (e.g., Vogelsberger et al. 2014; Schaye et al. 2015; Wang et al. 2015; Pillepich et al. 2018; Kannan et al. 2023, see Vogelsberger et al. 2020 for a review). A crucial factor in achieving these successes is the incorporation of sophisticated sub-grid models that approximately capture the complex baryonic physics involved, including stellar feedback originating within the multiphase ISM. However, these models face limitations in resolving the small-scale structure ( ≲ 1 kpc) of star-forming regions, as they do not explicitly model the physical processes within the ISM. While higher-resolution simulations enhance numerical accuracy, they struggle to naturally reproduce the observed complex multiphase ISM structure without incorporating more realistic feedback and ISM models (e.g. Agertz et al. 2013; Marinacci et al. 2019).  \nRecent advances in computational capabilities have enabled simulations of the cold ( ≲ 10 4 K) dense ( ≳ 1 cm -3 ) star-forming ISM in galaxies, attempting to resolve the physics including the cooling layers, H ii regions, Sedov-Taylor expansion, dust, and chemical processes, etc., with significantly enhanced fidelity to the underlying stellar feedback mechanisms (e.g. Hopkins et al. 2014, 2018, 2023; Rosdahl et al. 2015; Agertz & Kravtsov 2015; Marinacci et al. 2019; Kannan et al. 2020b). In these simulations, stellar feedback is typically modeled by averaging stellar properties over a stellar initial mass function (IMF), assuming that the IMF is fully sampled by stellar particles with masses of ≳ 10 3 M ⊙ . The feedback contribution of each stellar particle is then determined by the IMF-averaged properties rescaled by its mass. This IMF-averaged method is probably appropriate for massive galaxies, such as the Milky Way (MW). However, in dwarf galaxies, this approach tends to unrealistically 'smear' the stochastic mass, metal, and energy injections from the rare massive stars in space and time. This leads to an exacerbated 'over-cooling' problem, overly e ff ective photoionization feedback, and an incorrect baryonic cycle within these smaller galaxies (Su et al. 2018; Applebaum et al. 2020; Smith 2021).  \nFor the above reasons, some state-of-the-art simulations follow the so-called 'star-by-star' approach that explicitly samples individual massive stars from the IMF and tracks their feedback in simulations of both ISM boxes and entire dwarf galaxies (e.g., Revaz et al. 2016; Hu et al. 2017, H17; Emerick et al. 2019, E19;  \nApplebaum et al. 2020; Lahén et al. 2020; Smith et al. 2021; Hu et al. 2021; Gutcke et al. 2021; Rathjen et al. 2021; Calura et al. 2022; Andersson et al. 2023; Steinwandel et al. 2023; Lahén et al. 2023, L23). Notably, some of these dwarf galaxy simulations have reached unprecedented resolutions, albeit at the cost of modeling radiative feedback through approximate prescriptions. For example, Hu et al. (2016) simulated dwarf galaxies with a gas mass of 4 × 10 7 M ⊙ and a resolution of 4 M ⊙ by assuming a constant far-ultraviolet (FUV) radiation field. H17 improved upon their model by adopting optically thin prescriptions for a variable interstellar radiation field (ISRF) and Strömgren type approximation for photoionization feedback (Hopkins et al. 2012a, 2018) to study the nonlinear combination of feedback from radiation and SN explosions. Similar approximate treatments were then further developed and implemented by Hu (2019), Smith et al. (2021), and Steinwandel et al. (2023). Recently, based on the H17 feedback model, the GRIFFIN project simulated dwarf galaxies in isolation (Hislop et al. 2022; L23) and mergers (Lahén et al. 2019, 2020), providing systematic studies of various numerical and physical e ff ects on the formation and population of star clusters.  \nAlthough these models have demonstrated some successes in reproducing self-regulated star formation, multiphase ISM structures, galactic outflows, and cluster mass functions, their treatments of stellar feedback and ISM thermochemistry remain relatively simplified. Firstly, these models have calibrated e ff ective photoionization feedback with idealized tests of H ii region expansion (e.g. Bisbas et al. 2015). However, it is unclear whether these approximations for radiative feedback are su ffi ciently accurate in the wide variety of realistic environments encountered in simulations (Rosdahl et al. 2015; Fujimoto et al. 2019; Kannan et al. 2020a). The lack of long-range radiative e ff ects in these e ff ective models hampers their ability to accurately predict the thermal structure of CGM and IGM, leading to systematic errors in outflow properties and loading factors (e.g. Emerick et al. 2018). To more accurately capture stellar feedback and ISM thermochemistry, it is crucial to solve radiative transfer (RT) equations on the fly. Sophisticated numerical algorithms have been developed to make fully coupled radiation-hydrodynamics computationally feasible (e.g. Aubert & Teyssier 2008; Petkova & Springel 2009; Wise & Abel 2011; Rosdahl et al. 2013; Skinner & Ostriker 2013; Jaura et al. 2018; Kannan et al. 2019; Smith et al. 2020; Chan et al. 2021; Peter et al. 2023). A few models have already incorporated explicit RT solvers that, so far, are only a ff ordably realized for simulations of small dwarf galaxies with a modest resolution ( ≳ 100 M ⊙ ) or limited time coverage ( ≲ 500 Myr, e.g., E19; Agertz et al. 2020; Katz et al. 2022; Martin-Alvarez et al. 2023; Andersson et al. 2024). Additionally, these works may be limited by inadequate spatial and temporal resolution, which can result in inaccurate ionization feedback (Deng et al. 2024). Due to the improved e ffi ciency and accuracy of arepo -rt (Kannan et al. 2019; Deng et al. 2024), now on-thefly RT simulations are already obtainable for intermediate-sized dwarf galaxies (e.g., M gas = 4 × 10 7 M ⊙ ) with a high resolution of ∼ 1 M ⊙ .  \nIn addition, only a few models (e.g., E19) have taken into account the critical aspect of metallicity dependence on stellar properties such as lifetimes, photon production rates, mass-loss rates, and wind velocities. Notably, the ionizing photon production rate of OB stars tends to decrease with increasing metallicity (e.g., Schaerer 2002, S02; Lanz & Hubeny 2003, LH03), as metal-poor stars are generally more compact and hotter. Conversely, wind mass-loss rates are known to increase with metallicity (e.g. Vink et al. 2001; Vink & Sander 2021), a consequence  \nof boosted radiation pressure due to increased opacity in their atmospheres. These metallicity-dependent factors can introduce orders-of-magnitude di ff erences in the feedback budget of stars formed in either enriched or pristine environments, which must be reflected in cosmological simulations.  \nFinally, increases in simulation resolution must also be accompanied by more detailed cooling models capable of accurately capturing the multiphase structure of the ISM at smaller scales. The coupling between galactic star formation and outflow properties is highly sensitive to the ISM phase structure and the associated cooling rates, which have a complex dependence on metallicity, ultraviolet (UV) radiation fields, and nonequilibrium chemistry (e.g. Richings & Schaye 2016). However, Kim et al. (2023a) emphasized that the tabulated cooling models or fitting functions widely used for low-temperature ( < 10 4 K) cooling in many low-resolution galaxy formation simulations (e.g. Rosen & Bregman 1995; Smith et al. 2017; Hopkins et al. 2018; Ploeckinger & Schaye 2020) have significant discrepancies with those developed with detailed atomic or molecular physics by the ISM community (e.g. Wolfire et al. 2003; Gong et al. 2017; Bialy & Sternberg 2019; Kim et al. 2023a). Although the former models are computationally e ffi cient and capable of reproducing the cosmic star formation history and galaxy populations, they are inadequate at capturing the thermodynamics of the cold neutral ISM, especially when also modeling strong local radiation fields.  \nIn this paper, we introduce the Realistic ISM modeling in Galaxy Evolution and Lifecycles (RIGEL) model, a selfconsistent model for simulating stellar feedback and the multiphase ISM in dwarf galaxies. Our model brings together on-thefly RT, detailed photochemistry and cooling mechanisms, and metallicity-dependent feedback from individual massive stars. In this model, gravity, hydrodynamics, and the radiation field are explicitly modeled using the radiation hydrodynamics (RHD) code arepo -rt . The photochemistry model integrates nonequilibrium thermochemistry for H and He, along with the equilibrium abundances of key species like C i , C ii , O i , and CO, within irradiated ISM contexts. The physically motivated cooling model includes accurate prescriptions for all ISM phases, from hot ionized gas to cold molecular gas. Properties of massive stars, such as their lifetimes, photon production rates, mass-loss rates, and wind velocities, are determined by their initial masses and metallicities based on a library that incorporates a variety of stellar models. The remainder of the paper is organized as follows. In Section 2, we give an overview of the RIGEL model. In Section 3, we present the results of a suite of simulations of an isolated dwarf galaxy, exploring variations in numerical resolution (1 M ⊙ and 10 M ⊙ ), metallicity (0 . 02 Z ⊙ and 0 . 2 Z ⊙ ), and the e ff ects of turning RT on or o ff . We discuss the implications and caveats of our model in Section 4. Lastly, we give a summary in Section 5.", '2.1. Gravity, hydrodynamics, and radiative transfer': "The simulations presented in this work were performed with the moving-mesh hydrodynamic code arepo (Springel 2010; Pakmor et al. 2016, see Weinberger et al. 2020 for the public release 1 ). Gravity was solved with the tree particle-mesh method (Xu 1995; Bode et al. 2000; Bagla 2002), whereby the gravitational potential is split into short- and long-range contributions. The short-range gravity was calculated with a hierarchical  \noctree method (Barnes & Hut 1986), and the long-range gravity with the fast-Fourier-transformation method on a Cartesian mesh. Hydrodynamics was handled using a quasi-Lagrangian finite-volume method, which employs unstructured meshes generated via a Voronoi tessellation from discrete mesh-generating points. These volume-filling cells drift as they follow the local gas velocity, ensuring a pseudo-Lagrangian solution to the fluid dynamics. We employed a mass refinement scheme that maintains the cell masses around a target gas mass ( m gas) throughout the simulation domain by refining a cell into two if its mass exceeds twice the target mass or merging a cell with its neighbors if its mass falls below half of the target mass.  \nThe radiation field was explicitly modeled and coupled to gas hydrodynamics via the moment-based RHD solver arepo -rt (Kannan et al. 2019). The RT equations were solved by combining its zeroth and first moments with the M1 closure relation (Levermore 1984). The radiation spectra were discretized into seven broad bands, including infrared (IR, 0 . 1 -1 eV), optical (Opt., 1 -5 . 8 eV), far-ultraviolet (FUV, 5 . 8 -11 . 2 eV), Lyman-Warner (LW, 11 . 2 -13 . 6 eV), hydrogen ionizing (EUV1, 13 . 6 -24 . 6 eV), He i ionizing (EUV2, 24 . 6 -54 . 4 eV), and He ii ionizing (EUV3, 54 . 4 - ∞ eV) bands, to provide a comprehensive account for the photoionization or dissociation of various species and the dust-radiation coupling (see Table 1 for a summary and Section 2.4.2 for details). To avoid extremely small time steps, we use the reduced speed of light approximation (Gnedin & Abel 2001), where the selection of a reduced speed of light, ˜ c , is problem-dependent (e.g. Rosdahl et al. 2013). We also employed the spatial resolution correction methods outlined by Deng et al. (2024) to obtain convergent ionization feedback from the unresolved H ii regions. The coupled system of gravity, hydrodynamics, and RT was integrated using a two-stage secondorder Runge-Kutta scheme (i.e., Heun's method, Pakmor et al. 2016). Although omitted in this paper, magnetohydrodynamics implemented by arepo (Pakmor et al. 2011) can also be included in the simulations with RIGEL.", '2.2. Chemistry and cooling model': "Our model for chemistry and cooling closely follows Kannan et al. (2020b), with an update to the low-temperature ( ≲ 10 4 K) cooling by explicitly modeling the abundance of C and O species to better capture the thermodynamics in the warm and cold ISM. We tracked the time-dependent evolution of five primordial species - H, H + , H2, He, He + , and He 2 + -through arepo -rt . We refer the readers to Section 2.1 in Kannan et al. (2020b) for details of the primordial thermochemical network. When solving the nonequilibrium equations, we used a semi-implicit time integration approach based on the method outlined in Petkova & Springel (2009) and further improved by applying a photon absorption limiter to correct the temporal resolution issues due to large time steps (Deng et al. 2024). If the internal energy or one of the abundances changes by more than 10% during a time step, the equations will be solved implicitly by calling the functions from the SUNDIALS CVODE package (Hindmarsh et al. 2005).  \nWe also followed the enrichment and evolution of seven individual metal species including C, N, O, Ne, Mg, Si, and Fe. 2 The metals released by the stars di ff used and mixed numeri-  \ncally within the ISM, following the gas motion and mass exchange among gas cells as passive tracers. We defined the (absolute) metallicity, Z , of the ISM as the total mass fraction of these seven elements. Unlike Kannan et al. (2020b), which treats metal cooling as tabulated cooling rates from cloudy , we modeled high-temperature ( ≳ 10 5 K) and low-temperature ( ≲ 10 4 K) metal cooling separately. High-temperature metal cooling continued to be modeled using a pre-calculated cloudy table (see Appendix A.2). For low-temperature metal cooling, we tracked the equilibrium abundances of C, C + , CO, O, and O + to model the important metal cooling processes including fine structure lines of C ii 158 µ m, C i ∗ 610 µ m, C i ∗ 370 µ m, O i ∗ 63 µ m, and O i ∗ 146 µ m, as well as the CO rotational lines. The total abundances of these C and O species were directly obtained from the C and O abundances tracked by the gas particles. For dust involved in extinction and chemical reactions, we adopted a constant dust-to-metal ratio of 0.5, and thus Z g = Z d = 0 . 5 Z , where Z g and Z d are the gas and dust phase metallicity, respectively. Cosmic ray (CR) chemistry was also included by assuming constant ionization and heating rates (equation A.10) with 0.01 times the MW value. Specifically, the CR ionization rate for the neutral hydrogen atom is ζ cr,H i = 1 . 78 × 10 -18 s -1 , and ζ cr is used for ζ cr, H i if not specified.  \nTo model the cooling by the atomic and molecular lines of C and O species in the warm and cold ISM, we derived the equilibrium abundances of C, C + , CO, O, and O + , assuming that they are in a formation-destruction balance under the given abundance of primordial species and radiation field tracked by the RT module following the methods outlined by Kim et al. (2023a). For C, we calculated the equilibrium C + fraction by  \nx C + x C,tot = ζ pi,C + ζ cr,C ζ pi,C + ζ cr,C + α gr,C + n + α rr + dr,C + ne + k C + -H2 n H2 , (1)  \nwhere we have assumed that x C = x C,tot -x C + , since x CO is negligible in C / C + transition regions.  \nIn our model, C is ionized by photoionization and CRs. The photoionization rate, ζ pi,C = ζ Draine pi,C χ LW f s,C + 520 × 2 x H2 ζ cr, where χ LW is the LW band radiation intensity (tracked by arepo -rt ) normalized by the Draine (1978) value of 1 . 3 × 10 -14 erg cm -3 , ζ Draine pi,C = 3 . 43 × 10 -10 s -1 from Draine (1978), and f s,C is the self- and cross-shielding factor of C obtained following Tielens & Hollenbach (1985) and Gong et al. (2017). The second term accounts for the ionization by CR-generated LW-band photons in molecular gas, assuming that these photons are absorbed locally by gas (Heays et al. 2017). The CR ionization rate of C is proportional to that of H by ζ cr,C = 3 . 85 ζ cr. The rate coe ffi cients for the formation channels of C, including gas phase and grain surface recombination ( α rr + dr,C + and α gr,C + ) and the reaction of C + + H2 → CH + 2 ( k C + -H2 ), were taken from Gong et al. (2017).  \nThe O + abundance was determined by assuming x O + / x O,tot = x H ii because of the almost identical ionization potentials of O (13.62 eV) and H (13.60 eV) and the charge exchange between oxygen and hydrogen (see Section 14.7.1 in Draine 2011). Once the x C + and x O + were determined, we determined the CO fraction x CO by  \nx CO x C,tot = 2 x H2 GLYPH<2> 1 -max GLYPH<0> x C + / x C,tot , x O + / x O,tot GLYPH<1>GLYPH<3> 1 + ( n CO,crit / n ) 2 , (2)  \nabundances of the last seven elements, and we adopt the solar helium mass fraction Y = 0 . 2710.  \nwhere n CO,crit is the 'critical density' above which x CO / x C,tot > 0 . 5. We adopted the equilibrium fit in Gong et al. (2017),  \nn CO,crit = (4 × 10 3 Z ' d ζ -2 cr , 16 ) χ 1 / 3 LW       50 ζ cr , 16 Z ' 1 . 4 d       cm -3 , (3)  \nwhere ζ cr , 16 = ζ cr / (10 -16 s -1 ) and Z ' d is the dust metallicity normalized by the solar value.  \nBased on the chemical abundance of the main ISM coolants, we can calculate the cooling rate of di ff erent cooling channels. In our model, the cooling function is split into eight separate terms: primordial cooling from hydrogen and helium species ( Λ pri) tracked by the nonequilibrium chemical network and RT (local and ultraviolet background (UVB) photoionization heating is also incorporated into this term), high-temperature ( ≳ 10 5 K) metal collisional excitation lines cooling tracked by the cloudy look-up table ( Λ Z,CIE), nebular lines cooling in photoionized gas ( Λ Z,neb), equilibrium metal cooling in warm and cold gas ( Λ Z,C / O, including the cooling by C i , C ii , O i , and CO lines), cooling due to dust-gas-radiation interacting ( Λ dust), and heating due to photoelectric heating ( Γ PE) and CRs ( Γ cr). The net cooling rate ( Λ net) is then given by  \nΛ net = Λ pri( nj , N i γ , T ) + Z Z ⊙ Λ Z,CIE( T , ρ ) + Λ Z,C / O( nj , n C / O , N LW γ , T ) + Λ Z,neb( n H ii , ne , T , Z g) + Λ dust( ρ, T , D , N IR γ , Z d) -Γ PE( D , T , N FUV γ , Z g) -Γ cr( nj ) , (4)  \nwhere nj is the number density of all the primordial species ( j ∈ { H2, H i , H ii , He i , He ii , He iii } ), N i γ is the photon number density of all the photon bins from IR to He ii ionizing bands, T is the gas temperature, ρ is the density of the gas cell, and Z , Z g, and Z d are the total, gas-phase, and dust-phase metallicity, respectively. In Appendix A, we provide a detailed description of these cooling channels.", '2.3. Star formation and initial mass function sampling': 'In dwarf galaxies where the majority of stellar feedback comes from rare, massive stars, the stochasticity of when and where these massive stars form significantly influences the evolution of the galaxy. To capture this stochasticity, we developed a method to sample individual massive stars from a given IMF explicitly.  \nClosely following the method in Li et al. (2019), we identified gas particles as star-forming cells if they were cold ( T < T th), dense ( n H > n th), contracting ( ∇ · v < 0), self-gravitating ( |∇ · v | 2 + |∇ × v | 2 < 2 G ρ , Hopkins et al. 2013), and marginally Jeans-resolved ( m gas < f 3 J,n M J, where f J,n is a free parameter, M J is the local Jeans mass of the gas particle, see below). The choice of f J,n has large freedom, and typically this criterion is overwhelmed by the density and temperature threshold in our simulations. A star-forming gas cell can be converted to a star particle with a probability of P SF = ∆ t /τ ff , where τ ff = (3 π/ 32 G ρ ) 1 / 2 is the free-fall time of each cell and ∆ t is the length of the time step. We note that we assume the SFE per free-fall time, ϵ ff , equals unity as the resolution of our simulation is high enough to resolve the fragmentation of dense star-forming cores. For gas cells with m gas > f 3 J,s M J, we enforced instantaneous star formation to avoid the artificial fragmentation (Truelove et al. 1997), where f J,s is a free parameter similar to f J,n. We noticed there are di ff erent definitions of M J and selections of f J,s in the literature. We followed the numerical experiments by Grudi\'c et al. (2021),  \nwhich suggest that f J,s = 0 . 5 is acceptable with the Jeans mass defined as M J = π 1 . 5 c 3 s G -1 . 5 ρ -0 . 5 . The parameters used for the simulations in this work are presented in Section 3.1.  \nNewly created star particles inherit the mass and velocity of their progenitor gas cells ( m ⋆ = m gas,prog., v ⋆ = v gas,prog.) and are distinguished into two classes by a variable M fb. This variable determines whether the star particle actually hosts a M > M min,SN massive star and, if so, what its mass is. This feedback mass, M fb, typically di ff ers from the dynamical mass of the star particles, m ⋆ , and we determined the value of M fb using the method outlined below.  \nGiven the IMF, Φ ( M , Z ) ≡ d N / d M , and the minimum mass of massive stars (which can lead to core-collapse SN), M min,SN, we first drew a random number to determine whether the star particle hosts a massive star; each star particle with mass m ⋆ has a possibility, P host, of hosting a massive star,  \nP host = m ⋆ Z M max,IMF M min,SN Φ ( M , Z )d M / Z M max,IMF M min,IMF M Φ ( M , Z )d M , (5)  \nwhere M min,IMF and M max,IMF are the lower and upper mass limits of the IMF, respectively. If the star particle does not host an M > M min, SN star, we set M fb = -1. We note that in our highresolution simulations with m ⋆ ≲ 10 M ⊙ , P host < 1 always holds for the considered range of IMF variations.  \nIf a star particle does host a massive star, we drew another random number to determine its zero-age main-sequence (ZAMS) mass recorded in M fb from the inverse function of the cumulative function of IMF by M fb = F -1 ( y ), where y is a random variable with the value of [0 , 1] and F -1 is the inverse function of  \nF ( x ) ≡ Z x M min, SN Φ ( M , Z ) dM / "Z M max, IMF M min, SN Φ ( M , Z ) dM # . (6)  \nIf the IMF has a constant slope α (e.g., α = -2 . 3) at the M > M min,SN range, F -1 ( y ) has a simple form of  \nM fb = h ( M α + 1 max -M α + 1 min ) y + M α + 1 min i 1 / ( α + 1) . (7)  \nFor the simulations presented in this work, we adopted the IMF given by Chabrier (2003) and M min,IMF = 0 . 1 M ⊙ , M max,IMF = 100 M ⊙ , and M min,SN = 8 M ⊙ , respectively. We highlight here that implementing a variable IMF model (e.g., metallicitydependent IMF, see Kroupa et al. 2013 for a review) is straightforward under this framework and will be explored in future works.  \nWe emphasize that, in our high-resolution simulations, M fb usually exceeds m ⋆ when the star particle hosts a massive star. For example, a star particle with m ⋆ = 1 M ⊙ can host a 30 M ⊙ massive star that releases 23 M ⊙ ejecta when it explodes as an SN. To ensure mass conservation and avoid negative mass values of star particles, the star particles lose their own masses and release masses to the ambient medium in two distinct ways. Regarding the mass loss, we extracted the lost mass from all star particles continuously according to the IMF-averaged mass-loss rate through both the AGB winds and the SN explosion following the treatment outlined by Vogelsberger et al. (2013), regardless of whether it hosts a massive star. For mass release, star particles that do not host massive stars ( M fb = -1) can only return mass and metals to ISM mass through the AGB winds by 0 . 1 -8 M ⊙ stars. Only star particles with M fb ≥ 8 have the capability to release mass, metals, and energy through the SN explosion of the massive star they host and provide pre-SN feedback through ionizing radiation and stellar winds. As the actual mass  \nof star particles can be much smaller than the mass of ejecta, these particles do not consistently lose mass when they return mass by the stellar winds and SN ejecta. This inconsistency is counterbalanced by the mass loss resulting from the evolution of IMF-averaged stellar populations in all other star particles, ensuring mass conservation at the scale of star clusters.', '2.4. Zero-age main-sequence mass- and metallicity-dependent stellar feedback': "Once the massive stars are formed, we need to determine their lifetimes, photon production rates, mass-loss rates, and wind velocities based on their initial masses and metallicities. To build our library of these stellar properties, we gathered data on stellar evolutionary tracks, luminosities, mass-loss rates, and wind velocities from various references, as is detailed below, to encompass a range of metallicities from zero to the solar value and an initial mass range of M ∈ [8 -300] M ⊙ . We used the fitting formulae for ZAMS stellar properties in Tanikawa et al. (2020, hereafter T20) to bridge the properties in di ff erent literature and convert them into functions of M and Z . We found that the metallicity dependence of the ZAMS stellar radius and e ff ective temperature can be well approximated with a power law in the mass range of 8 -120 M ⊙ for Z ∼ 10 -6 -0 . 01 Z ⊙ , and that the luminosity is almost independent of Z , consistent with Larkin et al. (2023). Therefore, we did extrapolation with the power-law metallicity dependence to obtain the ZAMS stellar properties at Z = Z ⊙ , making a new reference model called T20 (EXT).  \nFor the convenience of implementation, we used a polynomial formula to fit the stellar properties, A ( M , Z ' ), where A can be the photon production rate in [s -1 ], mass-loss rate in [M ⊙ yr -1 ], wind velocity in [km s -1 ], and main-sequence time in [yr], as functions of the ZAMS mass at a given metallicity in the solar units Z ' ≡ Z / Z ⊙ and interpolated linearly in log 10 ( Z ' ) space. In practice, the interpolation over log 10 ( Z ' ) starts from a nonzero small metallicity Z 0, and the Z = Z 0 model is applied to all cases with Z < Z 0. We adopted Z 0 = 10 -8 Z ⊙ , below which stellar evolution is e ff ectively metal-free. 3 The polynomial fitting formula has the form  \nlog 10 A =        b 0 + b 1 x + b 2 x 2 + b 3 x 3 + b 4 x 4 , x ≤ x 0 , A 0 + d log 10 A dx GLYPH<12> GLYPH<12> GLYPH<12> x = x 0 ( x -x 0) , x > x 0 , (8) A 0 = b 0 + b 1 x 0 + b 2 x 2 0 + b 3 x 3 0 + b 4 x 4 0 , d log 10 A dx GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> x = x 0 = ( b 1 + 2 b 2 x 0 + 3 b 3 x 2 0 + 4 b 4 x 3 0 )  \nfor x ≡ log 10 ( M / M ⊙ ), where x 0 is a parameter chosen to optimize the fit quality at the high-mass end (which is better described by a linear relation between log 10 A and log M ). The relevant fit parameters, bi ( i = 0, 1, 2, 3, 4) and x 0, are shown in Appendix B (see Table B.1, Table B.2, and Table B.3). In most cases, we fixed b 4 = 0 to avoid overfitting.", '2.4.1. Stellar evolution': "In our model, we assume that the stellar properties do not evolve during their main-sequence lifetime. The main-sequence lifetimes of stars are determined by the metallicity-dependent function of Pop. II stars from Portinari et al. (1998) over Z =  \n0 . 0004 -0 . 05 combined with the function of metal-free Pop. III stars ( Z ≤ 10 -8 Z ⊙ ) from S02. In the left panel of Fig. 1, we present the main-sequence lifetimes, τ⋆ , as a function of the mass and metallicity of the main-sequence of zero age.  \nDuring the main-sequence stage, massive stars radiate photons (see Section 2.4.2 for details) and release mass and metals through stellar winds (Section 2.4.3). We assume that massive stars release unenriched gas by main-sequence winds and the yields by post-main-sequence winds and SN explosions are released instantaneously when the stars end their MS lifetimes.  \nOnce stars complete their main-sequence lifetime, they can release mass and metals through three channels: asymptotic giant branch (AGB) winds, core-collapse SNe, and type Ia SNe (SN Ia). Our AGB yield table covers the mass range of 1 -7 . 5 M ⊙ over Z = 0 . 0001 -0 . 02. This table was mainly adopted from Karakas (2010) and is supplemented by Doherty et al. (2014) and Fishlock et al. (2014) (see Pillepich et al. 2018 for details). The yield table for massive stars was adopted from the set M of Limongi & Chie ffi (2018, LC18), which covers the mass range of 13 -120 M ⊙ . These tables include the final integrated yields and the yields present in the stellar winds. The original LC18 tables have four metallicity bins expressed as [Fe / H] = -3 to 0. We converted the [Fe / H] values to Z = { 3 . 11 × 10 -5 , 3 . 11 × 10 -4 , 3 . 11 × 10 -3 , 0 . 0128 } based on the initial composition of C, N, O, Ne, Mg, Si, and Fe provided in their tables. A small fraction of < 8 M ⊙ low- and intermediatemass stars will explode as type Ia SNe following a delay-time distribution (see Section 2.4.5). We adopted the 'W7' model from Nomoto et al. (1997) for the SN Ia yields.  \nIn the rest of this section, we shall introduce our stellar radiation, main-sequence wind, and SN models in detail. We discussed the caveats and limitations of our stellar feedback model in Section 4.4.", '2.4.2. Radiation': "We included radiative feedback due to photodissosiation, photoionization, and photoheating of H2, H i , He i , and He ii species and the momentum transfer from photons to gas due to single and multiple scattering (i.e., radiation pressure). In our simulations, the star particles hosting alive M fb > 8 M ⊙ massive stars are sources of local radiation. We assume that the photon production rate, Q , remains constant at its ZAMS value throughout the main-sequence lifetime and drops to zero immediately after the main-sequence ends.  \nThe stellar spectra were discretized into seven bins (IR, Opt., FUV, LW, EUV1, EUV2, and EUV3) corresponding to our radiation bins. The photon production rate, Qi , for each spectral band was determined using polynomial fitting functions derived from the stellar mass-metallicity grid outlined later in this section. The mean energy per photon, mean ionization cross-section, photoheating rate, and the radiation pressure for each of the bins were calculated using a template synthesis spectrum taken from Bruzual & Charlot (2003) (see again Table 1 in Kannan et al. 2020b for details). Similar to Kannan et al. (2020b), we injected all photons into the nearest gas cell of the star particle to increase the probability of resolving the initial Strömgren mass. At time step ∆ t , this gas cell receives Qi ∆ t band i photons with mean energy ei . The direction of the photon flux, F i γ , was set to be radially outward from the star particle with a magnitude of | F i γ | = ˜ ceiQi ∆ t / V cell, where V cell is the volume of the gas cell receiving the photons.  \nThe photoionization and photon heating rates of species ' j ' are given by  \nζ pi , j = ˜ c X i ¯ σ i j N i γ (9)  \nand  \nΓ pi , j = ˜ c X i N i γ ¯ σ i j h i j , (10)  \nwhere ¯ σ i j and h i j are the mean ionization cross-section and photoelectron energy of photons in bin ' i ' interacting with species ' j ' (Kannan et al. 2019, 2020b).  \nRadiation can also directly couple momentum to gas by the single scattering. We added a source term in the momentum conservation equation of hydrodynamics and it is given by  \n∂ρ v ∂ t GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> rad. = 1 c X i F i γ         X j nj ¯ σ i j pij + κ i ρ ei         , (11)  \nwhere pij is the mean photon momentum and κ i is the opacity due to dust (see Table 1 in Kannan et al. 2020b for the values we adopted). The momentum boost due to multiple scattering was also modelled self-consistently by the method described in Kannan et al. (2019).  \nIn general, the photon production rate, Qi , increases with decreasing metallicity, as more metal-poor stars are more compact and hotter. To capture this trend, we constructed a metallicity grid for each photon energy band of ionizing photons with select models that satisfy a quasi-monotonic evolution of the total production rate of ionizing photons ( > 13 . 6 eV, combining EUV1, EUV2, and EUV3) with metallicity. Similarly, we selected models with monotonic evolution with metallicity for the Opt., FUV, and LW bands. For the IR band, we ignored the contribution of stellar IR radiation for simplicity, as it is employed for IR dustgas coupling (Kannan et al. 2019, 2020b). For the six photon energy bands above 1 eV, the stellar evolution models 4 used to construct the metallicity grids are described as follows.  \n- 1. For the Opt. and FUV bands, since the relevant feedback channels are unimportant for stellar feedback in metal-poor gas clouds, we calculated their production rates with the ZAMS stellar properties at Z = 10 -8 , 10 -6 , 10 -5 , 10 -4 , 0.01, and 1 Z ⊙ in T20 (EXT), assuming black-body spectra for simplicity.\n- 2. For the LW band, the metallicity dependence is unclear in the literature. For simplicity, we considered three models at Z = 10 -8 , 10 -4 , and 0 . 007 Z ⊙ in which Q decreases monotonically with Z . We adopted the results for Z ≥ 0 . 007 Z ⊙ from E19 (see their Appendix B) based on the OSTAR2002 stellar atmosphere models (LH03) and the ZAMS stellar properties from PARSEC (Bressan et al. 2012; Tang et al. 2014), which show minor variations with metallicity at Z ∼ 0 . 007 -1 . 2 Z ⊙ . At lower metallicities, we derived the photon production rates from the Z = 10 -4 Z ⊙ ZAMS properties in  \nLifetime  \nRadiation  \nWind  \nFig. 1. Stellar properties as a function of ZAMS mass and metallicity. Left: The main-sequence lifetime ( τ⋆ ) as a function of ZAMS mass and metallicity. The color of the curves from light to dark indicates the metallicity from low (10 -8 Z ⊙ ) to high (1 Z ⊙ ). Middle: Photon production rate ( Q ) of radiation bins from optical to He ii ionizing band as a function of ZAMS mass and metallicity. Similarly, the color of the curves from light to dark indicates the metal abundance from low (10 -8 Z ⊙ ) to high (1 Z ⊙ ). We shifted the curves of the Opt., FUV, LW, and EUV1 bands by multiplying constants of 10 4 , 10 3 , 10 2 , and 10, respectively, to avoid overlapping di ff erent bins. Right: Similar to left but for the wind mass-loss rate ( ˙ m w) and wind velocity ( v w). The color of the curves from light to dark indicates the metal abundance from low (10 -6 Z ⊙ ) to high (1 Z ⊙ ). The main-sequence stellar winds from extremely metal-poor stars with Z < 10 -6 Z ⊙ are negligible and not considered in our model.  \n<!-- image -->  \nWind  \nTable 1. Summary of the radiation energy bins adopted in our simulations.  \n- T20, assuming black-body spectra, and used the results for metal-free stars from Gessey-Jones et al. (2022, GJ22, see their Fig. 5) to cover Z ≤ 10 -8 Z ⊙ .\n- 3. For the EUV1 and EUV2 bands that dominate the ionization feedback from massive stars, we considered four models at Z = 10 -8 , 0 . 001 , 0 . 0282, and 1 Z ⊙ . We used the Z = 0 model in S02 (see their table. 6) for Z ≤ 10 -8 Z ⊙ . For Z = 0 . 001 Z ⊙ , we applied the OSTAR2002 grid (LH03) to the stellar parameters in T20 and took the average of the rates at ZAMS and the end of the main-sequence. Here, the stellar parameters at Z = 0 . 001 Z ⊙ were calculated by interpolating between the 10 -4 and 0 . 01 Z ⊙ models in T20. For Z = 0 . 0282 Z ⊙ , we used the fit in S02. For Z ≥ 1 Z ⊙ , we used the results from E19 at solar metallicity.\n- 4. For the EUV3 band, as optimistic estimates, we used the ZAMS stellar properties from T20 (EXT) to derive Q with black-body spectra for Z = 10 -8 , 10 -6 , 10 -5 , 10 -4 , 0.01, and 1 Z ⊙ .  \nWe summarize the properties of the seven radiation bands in Table 1. In the middle panel of Fig. 1, we present the photon production rate, Qi , as a function of ZAMS mass and metallicity for the six radiation bands above 1 eV with metallicity dependence (excluding the IR band).", '2.4.3. Winds of massive stars': 'Wemodeled the main-sequence winds of massive stars by injecting the mass, ˙ m w ∆ t , and momentum, √ 2 L w ˙ m w ∆ t , fluxes from a  \nstar to its nearest 32 neighboring gas cells weighted by their solid angle opening to the star (see Li et al. 2019 for details), where ˙ m w is the wind mass-loss rate and L w = 0 . 5 ˙ m w v 2 is the wind luminosity or power given the (terminal) wind velocity, v w.  \nSince the massive stars sampled in our simulations are typically hotter than the bi-stability-jump (with e ff ective temperatures T e ff > 30 000 K), we used the fit formula of Eq. 3 in Vink & Sander (2021) to calculate v w given the stellar luminosity, L , e ff ective temperature, T e ff , and metallicity, Z . The wind velocity was then combined with L , T e ff , Z , the stellar mass, M , and escape velocity, v esc = √ 2 GM / R ⋆ (given the stellar radius R ⋆ ), to derive the mass-loss rate, ˙ m w, through Eq. 24 in Vink et al. (2001). Similar to the calculation of photon production rates in the Opt. and FUV bands, we used the ZAMS values of M , L , T e ff , and R ⋆ from T20 (EXT) to calculate v w and ˙ m w for Z = 10 -6 , 10 -5 , 10 -4 , 0.01, and 1 Z ⊙ . Here, we ignored the weak main-sequence stellar winds from extremely metal-poor stars with Z < 10 -6 Z ⊙ for simplicity. The results were fit with Eq. (8), and the fit parameters are presented in Appendix B.We present the wind mass-loss rate, ˙ m v, and wind velocity, v v, as a function of ZAMS mass and metallicity in the right panel of Fig. 1.', '2.4.4. Core-collapse supernovae': "We assume that all stars with initial masses larger than 8 M ⊙ explode as core-collapse SNe and release 10 51 erg of thermal energy into their neighboring gas cells at the end of their main-  \nquence lives. We notice that the SN explosion energy has been modeled as a mass-dependent function in a few works (e.g. Gutcke et al. 2021; Steinwandel & Goldberg 2023), while its dependence on metallicity is still unclear (see Section 4.4 for details).  \nSince we adopted a 'pure thermal' scheme to inject SN energy, the feedback from SN can be significantly underestimated if we cannot resolve the Sedov-Taylor (ST) phase during which the thermal energy converts to kinetic energy while the total energy is conserved (i.e., the over-cooling issue, Katz 1992). To maximize the possibility of resolving the ST phase, we injected all the explosion energy to the nearest neighboring cell of the SN while distributing the ejected mass and metals to 32 neighbors weighted by the solid angle opening to the SN.  \nKim & Ostriker (2015) show that the resolution is su ffi -ciently high to ensure the momentum convergence if m gas,SN < M shell / 27, where m gas,SN is the mass of the gas cell that received the SN energy and M shell is the total swept-up mass at the end of the energy conserving ST phase (shell formation) given by  \nM shell = 1680 M ⊙ E 10 51 erg ! 0 . 87 GLYPH<18> n cm -3 GLYPH<19> -0 . 26 . (12)  \nFollowing equation (12), we can estimate the maximum density where the ST phase can be safely resolved, and it is 8 × 10 6 cm -3 for our fiducial 1 M ⊙ resolution and 10 3 cm -3 for 10 M ⊙ resolution. We have verified that the ST phase is resolved in most cases when the gas mass resolution is below 10 M ⊙ (see Appendix C, also see Hu 2019; Steinwandel et al. 2020, for the convergence tests in other numerical models showing similar results).", '2.4.5. Type Ia supernovae': "For type Ia SNe, we adopted a delay-time distribution (DTD) formalism to compute the explosion probability of each star particle as a function of time. In the DTD formalism, the type Ia rate is determined by  \n˙ N Ia( t ) = N Ia Z t 0 ˙ M ⋆ ( t ' ) Ψ ( t -t ' )d t ' , (13)  \nwhere N Ia is an observable, ˙ M ⋆ ( t ) is the star formation rate (SFR), and Ψ ( t ) is the DTD. For each star particle in our simulations, ˙ M ⋆ ( t ) = m ⋆δ ( t -t 0), where δ ( t ) is the Dirac Delta function, and t 0 is the formation time of the star particle. For simplicity, we denote t 0 = 0. This gives a simplified SN Ia rate for each star particle: ˙ N Ia( t ) = m ⋆ N Ia Ψ ( t ).  \nIn our high-resolution simulations, each star particle should at most have one SN Ia. Similar to our IMF-sampling scheme, we first drew a random number to determine whether the star particle is an SN Ia progenitor. Each star particle with mass m ⋆ has a possibility, P Ia, of being an SN Ia progenitor,  \nP Ia = m ⋆ N Ia . (14)  \nIf this star particle is an SN Ia progenitor, we drew another random number, y , to determine its explosion time from the inverse function of the cumulative function of DTD by t = F -1 ( y ), where  \nF ( t ) = Z t 0 Ψ ( t ' )d t ' . (15)  \nWeadopted a power-law DTD model that consists of a power law in time,  \nΨ ( t ) = ( 0 , t < τ⋆ (8M ⊙ , Z ) At -s , t > τ⋆ (8M ⊙ , Z ) , (16)  \nArticle number, page 8 of 30  \nwhere A is a normalization factor such that the integration of Ψ ( t ) over a Hubble time, t H, equals 1 (i.e., R t H 0 Ψ ( t ' )d t ' = 1) and τ⋆ (8M ⊙ , Z ) is the lifetime of an 8 M ⊙ massive star with an initial metallicity of Z . Specifically, we adopted s = 1 . 12 and N Ia = 1 . 3 × 10 -3 SNM -1 ⊙ (Maoz et al. 2012; Vogelsberger et al. 2013). While our model accounts for it, we leave the discussion of SN Ia feedback and enrichment for our future cosmological simulations, as it occurs too late to influence the cloud-scale star formation that is the focus of this paper.", '3.1. Initial conditions and model parameters': "To evaluate the performance of the RIGEL model, we ran a suite of simulations of an isolated dwarf galaxy. We generated the isolated disk initial conditions using the M ake D isk code developed by Springel et al. (2005). We created a disk with parameters identical to the fiducial initial condition used in H17, which is also extensively used in recent isolated dwarf simulations (e.g., Gutcke et al. 2021; Hislop et al. 2022; L23). The dark matter halo follows a Hernquist profile with an NFW-equivalent concentration parameter, c = 10. The virial mass of the halo is M vir = 2 × 10 10 M ⊙ , the virial radius is R vir = 44 kpc, and the spin parameter is λ = 0 . 03. The rotation-supported galaxy disk has an initial stellar disk with a mass of 2 × 10 7 M ⊙ and a scale length of 0 . 73 kpc and an initial gas disk with a mass of 4 × 10 7 M ⊙ and a scale length of 1 . 46 kpc. The scale height of both the gas and stellar disks is 0.35 kpc. The star particles in the initial disk only interact with gas via gravity but do not contribute to stellar feedback or enrichment.  \nThe gravitational softening lengths for dark matter (DM) and star were determined following Hopkins et al. (2018) by ϵ DM = 1 . 7 pc m 1 / 3 DM ( ρ DM / 10 9 M ⊙ kpc -3 ) -1 / 3 and ϵ⋆ = 0 . 34 pc m 1 / 3 ⋆ ( ⟨ n SF ⟩ / 10 3 cm -3 ) -1 / 3 , where ⟨ n SF ⟩ is the mean density for star formation. We list the softening lengths used in this work in Table 2. The softening of gas cells was determined adaptively with a minimum length of 0 . 004 pc. We set the reduced speed of light as ˜ c = 1000 km s -1 to capture the propagation of the ionization front in a typical ISM environment with a reasonable time lapse (Rosdahl et al. 2013).  \nTo avoid initial artificial starbursts, we relaxed the initial conditions by injecting thermal energy following the method outlined in H17 for 500 Myr to create a turbulent ISM within the gas disc. In production simulations, we adopted a spatially uniform ionizing UV background radiation field from FaucherGiguère et al. (2009) with corrections for self-shielding in dense ISM (Rahmati et al. 2013). In addition, we set a background FUV field of G 0 , UVB = 0 . 00 324 for the C i ionization and photoelectric heating, which is the UVB value adopted from Haardt & Madau (2012) integrated over the 5.8-13.6 eV energy range. The star formation criteria used for production simulations are T th = 100 K, n th = 10 4 cm -3 for 1 M ⊙ resolution and n th = 10 3 cm -3 for 10 M ⊙ resolution. This choice of this density threshold was motivated by the observational implications (e.g. Parmentier et al. 2011; Lada et al. 2012). The necessary and sufficient Jeans numbers for star formation were set as f J,n = 0 . 1 and f J,s = 0 . 5.  \nIn this paper, we present a total of six simulations with different combinations of feedback channels, initial metallicities, and resolutions. We list the naming conventions and simulation setups in Table 2. We refer to the full RIGEL physics model as the 'full' model. In this paper, the full RIGEL physics simulation with 1 M ⊙ baryon resolution initialized with Z = 0 . 02 Z ⊙  \nis our fiducial run, named as '0 . 02Z ⊙ / full.' For comparison, we simulated the same galaxy without local radiative feedback (but with UVB), referred to as '0 . 02Z ⊙ / noRT.' We also simulated a galaxy initialized with 0 . 2 Z ⊙ with the 'full' physics, and it is referred to as the '0 . 2Z ⊙ / full' run. For the convergence test, we also ran each simulation with a lower resolution of 10 M ⊙ and these runs are referred to as 'low.'  \nWe ran all six simulations for 1 Gyr, which corresponds to roughly six orbit times at a 1 kpc galactocentric radius. As our simulated dwarf galaxies exhibit regular rotation and steady star formation activity, each ∼ 200 Myr time segment of the simulation can be considered as a realization of di ff erent random seeds. We additionally ran the 0 . 02Z ⊙ / low simulation using another random seed, which demonstrates minor variations in the global SFR and ISM structures. Since di ff erent disks spend different times to cool down and form the first star from the initial condition, we define the time origin of each simulation as the moment its first star forms.  \nFor each 1 M ⊙ -resolution simulation, ∼ 1 . 3 × 10 6 core-hours are required for every 1 Gyr of simulated time using the Intel Xeon Platinum 8160 ('Skylake') CPUs; the 10 M ⊙ -resolution simulations, which have ten times fewer particles and roughly 10 1 / 3 larger time steps, reduce the computational cost to ∼ 3 × 10 4 core-hours per 1 Gyr of simulation.  \nWe notice that though the disk properties of our galaxies are the same as those in H17 and L23, the initial metallicity of H17 is 0 . 15 Z ⊙ in our normalization, slightly lower than our 0 . 2 Z ⊙ galaxies, while the initial metallicity of L23 is 0 . 015 Z ⊙ , also slightly lower than our 0 . 02 Z ⊙ galaxies. Nonetheless, we think our results can be directly compared with them because such di ff erences are minor.  \nWe use cylindrical coordinates R and z to describe our simulations, where R is the galactocentric radius and z is along the rotation axis of the disc. We define the 'ISM region' as the R < 2 kpc and | z | < 1 kpc region.", '3.2. Star formation rates and the Kennicutt-Schmidt relation': "We first analyze the star formation properties of our simulations. The left panel in Fig. 2 shows the SFR as a function of simulation time, representing the star formation histories of the simulated galaxies, and the right panel shows the cumulative mass of newly formed stars in our simulations. In each simulation, there is a rapid rise in the SFR initially, followed by the SFR stabilizing at a self-regulated state due to the presence of feedback. Despite the average SFRs remaining relatively stable over long timescales, there are noticeable bursty fluctuations at shorter scales. In particular, the SFRs in the full feedback models exhibit fluctuations exceeding two orders of magnitude.  \nThe SFRs show an evident dependence on metallicity. The median SFR of 0 . 02Z ⊙ / full is only ∼ 6% of that in 0 . 2Z ⊙ / full. This phenomenon is mainly attributed to the low cooling rate, which suppresses the formation of cold, dense gas. Consequently, once the gas is heated to ≳ 10 4 K, a longer period is required for the metal-poor gas to cool down and convert to a star-forming state. Although the total feedback budget (which scales with the SFR) is relatively limited in low-metallicity environments, it becomes more di ffi cult for the ISM to revert to its cold phase once heated to its warm-di ff use phase due to ine ff ective cooling mechanisms and suppression of thermal instability (e.g. Schaye 2004; Walch et al. 2011; Bialy & Sternberg 2019).  \nTurning o ff radiative feedback has a larger impact on the SFR. The median SFR of 0 . 02Z ⊙ / noRT is 92 times higher than that of 0 . 02Z ⊙ / full and even five times higher than 0 . 2Z ⊙ / full.  \nWe will see in Section 3.4 that this di ff erence is due to the fact that 0 . 02Z ⊙ / noRT contains the highest amount of cold and molecular gas, as there is no ISRF present to heat the di ff use gas and keep it warm and atomic.  \nWe note that the SFRs exhibit good convergence with resolution. In all models, the cumulative stellar masses of the lowresolution simulations are slightly lower, but the di ff erence in stellar mass formed between the low- and high-resolution simulations is within a factor of two. The 0.02Z ⊙ / noRT simulations show the most significant di ff erence in cumulative stellar mass between the low- and high-resolution simulations, because the SFR of the high-resolution run is systematically higher than that of the low-resolution run before t ≈ 500 Myr. These systematic di ff erences can be attributed to some specific features of initial conditions for high- and low-resolution simulations, as we stochastically injected thermal energy to relax them (Section 3.1).  \nTo assess the performance of our model in regulating star formation, we analyze the Kennicutt-Schmidt (KS) relation (Kennicutt 1998) which establishes a connection between the surface densities of gas ( Σ gas) and the SFR ( Σ SFR) in our simulations. We also analyze the extended KS relation which substitutes the total gas surface density by including the contribution of stars' gravity ( Σ ⋆ , e.g. Shi et al. 2011, 2018). In Fig. 3, we present the KS relation ( Σ gasΣ SFR, left panel) and the extended KS relation ( Σ 0 . 5 ⋆ Σ gasΣ SFR, right panel) of the simulated highresolution galaxies for the 0 . 02Z ⊙ / full (red), 0 . 2Z ⊙ / full (green), 0 . 02Z ⊙ / noRT (blue) model. The SFR surface densities are measured as the average over 1 Myr within a series of 100 pc annuli from the simulation time t = 300 Myr to 1000 Myr. We also plot the mean Σ SFR of the [0,25], [25,50], [50,75], and [75,100] Σ gas ( Σ 0 . 5 Σ gas) percentile bins to show the trend.  \n⋆  \nOn the Σ gasΣ SFR plane (Fig. 3 left panel), the results from both full feedback simulations are in good agreement with the observed dwarf galaxies in the FIGGS survey measured with 400 pc resolution (Roychowdhury et al. 2015), whose 16-84 and 5-95 percentile ranges are presented as the dark blue and light blue shaded regions. For a specific Σ gas bin, the scatter in Σ SFR is 1 . 0-1 . 4 dex, which is comparable to FIGGS's scatter of 1 . 31 . 8 dex for the 0 . 15 < log( Σ gas) < 0 . 75 range. The medians of 0 . 2Z ⊙ / full data roughly follows the Elmegreen (2015) star formation relation for dwarf galaxies and outer disks, characterized by a power-law index of 2 (black dashed line). This power-2 star formation relation suggests that the thickness of the disk is regulated by the balance between self-gravity and vertical pressure, in regions with low stellar surface density (e.g. Ostriker et al. 2010; Elmegreen 2018). The 0 . 02Z ⊙ / full also roughly follows the power-2 star formation relation, while the intercept is smaller, implying a longer dynamical time-scale for star formation in the low metallicity environments due to the ine ffi cient cooling. The 0 . 02Z ⊙ / noRT results lie outside the 16-84 percentile range of the FIGGS data, exhibiting higher SFR surface density, but still lower than the classical Kennicutt (1998) relation. The slope of its KS relation is even shallower than the classic slope of 1.4, indicating a more rapid and e ffi cient star formation process in the absence of momentum injection and photoheating from radiative feedback.  \nOn the Σ 0 . 5 ⋆ Σ gasΣ SFR plane (Fig. 3 middle panel), 0 . 2Z ⊙ / full roughly follows the extended KS relation given by Shi et al. (2018). All three simulations exhibit a similar slope, indicating the significant role of pre-existing disk stars in regulating star formation in our simulated galaxies.  \nTable 2. Overview of the simulations reported in this work and simulation parameters.  \nNotes. The columns from left to right are the model name (first column), initial metallicity (second column), target mass of gas particles (third column), mass of star particles (both background and active stars, fourth column), mass of dark matter particles (fifth column), max softening length of gas (sixth column), softening length of stars (seventh column) and DM (eighth column), and the stellar feedback channels included in the models (ninth column). Abbreviations: PI: photoionization, RP: radiation pressure, PE: photoelectric heating, WD: stellar wind, SN: supernova.  \n<!-- image -->  \nFig. 2. Star formation histories of the simulated dwarf galaxies. Left: SFR as a function of time. The SFR is averaged over a time-scale of 10 Myr. The solid curves are the results from the high-resolution runs while the dotted curves are from corresponding low-resolution runs. Right: Cumulative mass of newly formed stars in our simulations. The SFRs exhibit good convergence with resolution. The 0.02Z ⊙ / noRT simulations show the most significant di ff erence in cumulative stellar mass between the low- and high-resolution simulations; nonetheless, the di ff erence remains within a factor of two.  \n<!-- image -->  \nThe results in the KS and extended KS planes of our full feedback simulations demonstrate the ability of the RIGEL model to capture the self-regulated galactic star formation environments. In the 0 . 2Z ⊙ / full galaxy, both the KS and extended KSrelations are consistent with the observed empirical relations, while the results of 0 . 02Z ⊙ / full suggest a longer dynamical timescale for star formation (a smaller intercept for the relations) in extremely metal-poor environments.", '3.3. Morphologies': 'To provide a panoramic view of the entire galactic disk, Fig. 4 presents the face-on surface density distribution (the first row) and mid-plane slices of gas temperature (the second row) and metallicity increment ( ∆ Z = Z -Z init, the third row) of the simulated galaxies at 900 Myr for the 0 . 02Z ⊙ / full, 0 . 2Z ⊙ / full, and 0 . 02Z ⊙ / noRT models. Young stars (less than 100 Myr old) are color-coded, ranging from purple for the youngest (0 Myr) to red for the oldest (100 Myr).  \nIn all three galaxies, the ISM exhibits a multiphase structure. The galactic disks display filamentary and clumpy patterns, with many bubbles amidst the clouds, which illustrate the ongoing stellar feedback processes. Bubbles with moderate density contrast are mainly carved out by early feedback from radiation and stellar winds, while those with high density contrast are primarily the result of SN explosions.  \nThe 0 . 02Z ⊙ / full galaxy presents a relatively smooth gas distribution and less clumpiness. Cold ( < 100 K) gas is only found in sporadic small clumps, while a large part of the disk remains warm. Hot gas also only appears in relatively isolated bubbles blown by the SN explosion. The smaller warm bubbles are expanding H ii regions of young massive stars. In contrast, in 0 . 2Z ⊙ / full, although warm gas still dominates, the structure is more filamentary and clumpy due to more rapid cooling caused by the higher metal content. Hot SN bubbles are interconnected to form superbubbles on a kiloparsec scale. The ISM in 0 . 02Z ⊙ / noRT is ripped into fragments by intense SN feedback due to its highest SFR (see Section 3.2). Interestingly, though a large fraction of the disk is occupied by hot superbubbles, a con-  \n)  \n]  \n2  \nc  \np  \nk  \n1  \nr  \ny  \nM  \n[  \n/  \nR  \nF  \nS  \n(  \n0  \n1  \ng  \no  \nl  \n2  \n3  \n4  \n5  \n6  \n7  \ngas-  \n1 Gyr  \n10 Gyr  \n100 Gyr  \n1000 Gyr  \n1  \n0  \n1  \n2  \nlog10(  \ngas / [M  \npc  \n2  \nFig. 6 quantifies the time evolution of the mass fractions (top row panels) and volume fractions (bottom row panels) of gas in the ISM region ( R < 2 kpc and | z | < 1 kpc) for the cold neutral medium (CNM, T < 100 K), the warm neutral medium (WNM, 100K < T < 8000 K), the warm ionized medium (WIM, 8000 K < T < 10 000 K), and hot ionized medium (HIM, T > 10 000 K) phases.  \n<!-- image -->  \nKennicutt 98  \nElmegreen 15  \nFIGGS  \nTHINGS  \nLITTLE THINGS  \n2  \n])  \nFig. 3. Star formation relations in the simulated dwarf galaxies. Left: Kennicutt-Schmidt relation for the simulated high-resolution galaxies for the 0 . 02Z ⊙ / full (red), 0 . 2Z ⊙ / full (green), and 0 . 02Z ⊙ / noRT (blue) model. The translucent data points are 1 Myr average SFR within a series of 100 pc annuli measured from the simulation time t = 300 Myr to 1000 Myr. The data with errorbars are the median Σ SFR of the [0,25), [25,50), [50,75), [75,100] gas surface density percentile bins, and the errorbars of Σ SFR show the 16 and 84 percentiles of the SFR surface density in each bin. The black solid line shows the standard KS relation from Kennicutt (1998); the black dashed line shows the power-2 relation for outer disk regions and dwarf irregular galaxies from Elmegreen (2015). The observational data collected by Elmegreen (2015), including the far outer regions of spirals and dwarf galaxies from the THINGS survey (Bigiel et al. 2010) and the dwarf galaxies from the LITTLE THINGS survey (Elmegreen & Hunter 2015), are shown as grey and black dots, respectively. The blue-shaded region is the 5-95 and 16-84 percentile range of the observed dwarf galaxies in the FIGGS survey (Roychowdhury et al. 2015). Right: extended KS relation ( Σ SFR ∝ ( Σ 0 . 5 ⋆ Σ gas) 1 . 09 , Shi et al. 2011, 2018) measured in the same way; the black solid line is the fitting of observed data of star-forming galaxies of various types from Shi et al. (2018).  \nsiderable amount of cold gas coexists with these large hot bubbles (see Section 3.4 for the phase diagrams and distributions of gas). This coexistence is explained by the fact that winds and SNe primarily heat and disrupt molecular clouds locally, while maintaining the warm state of di ff use gas requires FUV radiation. Furthermore, the absence of radiation prevents e ffi cient evacuation of the dense gas surrounding massive stars, which suppresses the e ffi ciency of SN feedback. Consequently, the cold molecular gas can persist instead of being destroyed by FUV radiation that permeates the entire disc.  \nThe metallicity distribution in these galaxies shows significant spatial variations, where areas that are enriched beyond normal levels are associated with SN bubbles. Small-scale eddies can be found that trace the turbulent flow of gas. Generally, the central regions exhibit higher metallicity due to the higher SFR, while metallicity decreases gradually as distance from the galactic center increases. In the non-star-forming disk, the metallicity distribution appears relatively uniform and smooth, suggesting that the over-enriched bubbles could blend with the ISM and fade away rapidly.  \nFig. 5 presents the same projections and slices but viewed from an edge-on perspective. In the 0 . 02Z ⊙ / full galaxy, the disk exhibits a spindle shape that flares outwards in the outer disk due to pressure support. Enriched gas is confined to specific chimneys within its CGM, while a significant portion of the CGM remains pristine. On the other hand, both the 0 . 2Z ⊙ / full and 0 . 02Z ⊙ / noRT galaxies feature a more enriched CGM, with the enriched gas occupying almost the entire volume of the plotted box. The spindle shape is disrupted by multiple feedback-driven chimneys through which hot outflow gas escapes. Notably, the  \nCGM of the 0 . 02Z ⊙ / noRT galaxy is hotter and more enriched than that of the 0 . 2Z ⊙ / full galaxy. These results demonstrate that the inclusion of energetic radiative feedback reduces the outflow of gas and metals by suppressing the SFR. Further analysis of the outflow properties of our galaxies will be detailed in Section 3.6.', '3.4.1. Mass and volume fractions of different interstellar medium phases': 'In all three simulations, warm gas dominates both the mass and volume fractions. The majority of the mass is found in the WNM phase, while WIM occupies a smaller mass but a similar volume in the ISM as the WNM. The CNM only makes up about 10 -3 of the ISM mass in the 0 . 02Z ⊙ / full galaxy, with its mass fraction increasing by a factor of around 10 in the 0 . 2Z ⊙ / full galaxy due to more e ff ective metal cooling. Notably, in the 0 . 02Z ⊙ / noRT galaxy, the CNM mass fraction is about 40 times higher than in the 0 . 02Z ⊙ / full galaxy in the absence of radiative feedback. As a result, 0 . 2Z ⊙ / full and 0 . 02Z ⊙ / noRT exhibit higher SFRs due to the availability of more cold, dense gas for star formation. The higher SN rates in these two galaxies pro-  \nSFR', '0.2Z0 /full': 'Fig. 4. Face-on images of the integrated gas mass surface density ( Σ gas; top row), along with mid-plane slices of gas temperature ( T ; middle row), and metallicity increment ( ∆ Z = Z -Z init; bottom row) of the simulated galaxies at 900 Myr. The three columns correspond to the di ff erent runs, as labeled above the panels. From left to right, they are 0 . 02Z ⊙ / full, 0 . 2Z ⊙ / full, and 0 . 02Z ⊙ / noRT, respectively. The overset points represent young stars ( < 100 Myr) with ages color-coded from violet (youngest) to red (oldest).  \n<!-- image -->  \nduce more hot gas as well. In 0 . 02Z ⊙ / full, the volume fraction of HIM is only ∼ 10 -3 due to the low SN rate, while it constitutes a significant volume fraction ( ∼ 0 . 1) in 0 . 2Z ⊙ / full and 0 . 02Z ⊙ / noRT. On the other hand, the volume occupied by CNM is always negligibly small in all three galaxies as CNM only appears in dense clumps. The mass and volume fractions of the low-resolution simulations are presented in dotted curves, which show very similar warm and cold gas factions to those of the high-resolution simulations. The HIM volume fractions exhibit the largest systematic di ff erences with resolution, possibly due to the slightly lower SFR and potentially unresolved cooling at  \nlow resolution, which results in values up to 0.8 dex lower in the low-resolution simulations.', '3.4.2. Phase diagram and interstellar medium distributions': 'To better understand the ISM structures in our simulated galaxies, we present the gas phase diagrams for the high-resolution galaxies averaged over the snapshots from t = 300 Myr to t = 1000 Myr with a time interval of 10 Myr in Fig. 7. We also show the median temperature (red dashed curves) and the equilibrium temperature (grey curves) as a function of density. The median curves roughly follow the equilibrium curves, though  \n0.2Z0 /fullFig. 5. Similar to Fig. 4, but viewed from an edge-on perspective.  \n<!-- image -->  \nthere is a significant dispersion in the gas distribution due to the feedback processes and the turbulence stirring. The gas distribution presents a typical bi-stable structure, featuring stable warm and cold phases connected. The hot gas heated by SN explosions lies in the upper left quadrant of the phase diagrams, and the gas in the diagonal branches of the lower right quadrant represents thermally unstable gas transitioning from its warm to cold phase. In the 0 . 02Z ⊙ / full galaxies, both the upper left and lower right quadrants are less populated due to the lower SN and cooling rates. The photoionized gas within H ii regions manifests as a faint, narrow, ∼ 10 4 K isothermal line in 0 . 02Z ⊙ / full and 0 . 2Z ⊙ / full galaxies.  \nIn the first column of Fig. 8, we present the mass-weighted probability distribution functions (PDFs) of gas density (top) and temperature (bottom) in the ISM region of the simulated  \ngalaxies. The dotted lines are the results of the low-resolution simulations, which closely resemble the high-resolution ones. Across all three galaxies, the predominant gas component is warm ( ∼ 10 4 K) and di ff use ( ∼ 1 cm -3 ). The distributions of warm gas in 0 . 02Z ⊙ / full and 0 . 2Z ⊙ / full exhibit substantial similarity, suggesting that the thermal structure of warm gas is insensitive to metallicity. This is due to the dominance of Ly α cooling over the metallicity-dependent cooling and heating processes in this phase. When the density exceeds ∼ 1 cm -3 , the gas in 0 . 2Z ⊙ / full is more susceptible to thermal instability, leading to easier cooling and condensation into cold, dense gas. Thus, 0 . 2Z ⊙ / full displays a distinct two-phase behavior with an additional peak of cold gas at temperatures ≲ 10 2 K. As a result, the 0 . 2Z ⊙ / full galaxy contains significantly more cold, dense  \nFig. 6. Time evolution of the mass fractions (top row panels) and volume fractions (bottom row panels) of gas in the ISM regions ( R < 2 kpc and | z | < 1 kpc) in the cold neutral medium (CNM, T < 100 K, blue curves), the warm neutral medium (WNM, 100 < T < 8000 K, green curves), the warm ionized medium (WIM, 8000 < T < 100 000 K, orange curves), and hot ionized medium (HIM, T > 100 000 K, red curves) phases. The dotted curves are from the low-resolution runs. Fig. 8 shows how these phases are distributed across the mass-weighted ISM probability distribution functions. Both higher metallicity and noRT runs exhibit significantly more CNM and HIM fractions due to the more e ffi cient gas cooling and collapse (induced by higher cooling rate and lower heating rate in 0 . 2Z ⊙ / full and 0 . 02Z ⊙ / noRT, respectively), and the resultant higher SFRs.  \n<!-- image -->  \ngas than 0 . 02Z ⊙ / full. Additionally, the 0 . 2Z ⊙ / full also contains more hot, low-density gas because of its higher SN rate.  \nNotably, the ISM phase structure of the 0 . 02Z ⊙ / noRT galaxy di ff ers significantly from that of 0 . 02Z ⊙ / full, even for the warm di ff use gas. As we have noticed from the morphologies, large superbubbles and cold, dense clumps are coexistent in this galaxy. This feature is also reflected in the PDFs, where the fractions of both hot and cold gas are notably elevated compared to the other two galaxies. This discrepancy arises because the absence of heating by the ISRF allows the di ff use ISM to cool more easily, rendering it susceptible to thermal instability. As a result, it contains a higher proportion of cold dense gas compared to 0 . 02Z ⊙ / full, despite its lower metallicity and higher SFR. Furthermore, 0 . 02Z ⊙ / noRT exhibits the highest amount of hot gas. This is not only because its higher SFR leads to more SN explosions, but also because the more clustered star formation in the absence of radiative feedback increases the clustering of SNe (see Section 3.5 and Smith et al. 2021).  \nNow, we shift our focus to the distribution of hydrogen species tracked by the nonequilibrium chemical network across its ionized, neutral, and molecular phases (H ii , H i , and H2). In the second to fourth columns of Fig. 8, we present the mass distribution of gas for di ff erent galaxies weighted by the H ii (purple), H i (orange), and H2 (blue) fractions. We see that all three galaxies are atomic-dominated. In these galaxies, H i gas dominates at densities ranging from n H ∼ 10 -2 . 5 cm -3 to at least 10 3 cm -3 , covering a broad range of temperatures of T ∼  \n10 -10 4 K. H i can be ionized to H ii by collisional ionization in warm and hot gas with T ≳ 20 000 K, by UVB in unshielded low-density regions ( ≲ 10 -2 cm -3 ), and by local photoionization in the vicinity of young massive stars. H ii gas dominates the warm di ff use ISM and the hot low-density gas that corresponds to SN-driven bubbles, and it has an extended distribution towards the dense ( ∼ 10 cm 3 ) cold ( ∼ 100 K) gas. In the 0 . 02Z ⊙ / noRT galaxy, the distribution of H ii gas exhibits symmetry in density in the absence of radiation fields, whereas in the full feedback simulations, particularly in 0 . 2Z ⊙ / full, there is a departure from symmetry with a high-density tail comprising the photoionized gas within H ii regions.  \nThe H2 fractions in the full feedback galaxies are primarily determined by competition between dust-phase formation ( ∝ Z d) and photodissociation by LW radiation. On the other hand, in the noRT galaxy, H2 can only be destroyed slowly by collisional processes and CR ionization. In the 0 . 02Z ⊙ / full galaxy, H2 remains a minority component even in dense gas with > 10 3 cm 3 due to its low metallicity. Half of the total H2 mass exists in n H < 1 cm -3 di ff use gas where the H2 fraction is low. The situation is similar in the 0 . 2Z ⊙ / full galaxy. Although there is more cold, dense gas where H2 can form, the H i -H2 transition occurs at high densities of > 10 3 cm -3 with a significant fraction of the H2 mass existing in H i -dominated gas. The 0 . 02Z ⊙ / noRT galaxy exhibits a broader distribution of H2 across large temperature and density ranges. H2 even exists in warm-hot gas ( T ∼ 10 4 -10 6 K) with densities n H ∼ 10 -4 -10 -2 cm -3 , because the gas-phase for-  \n0  \nFig. 7. Mass-weighted gas phase diagrams for the high-resolution galaxies averaged over the snapshots from t = 300 Myr to t = 1000 Myr with a time interval of 10 Myr. The red dashed curves represent the median temperature within each density bin. The grey curves are the equilibrium temperature assuming G 0 = 0 . 01, 0 . 1, and 0 . 00324 for 0 . 02Z ⊙ / full, 0 . 2Z ⊙ / full, and 0 . 02Z ⊙ / noRT, respectively.  \n<!-- image -->  \nFig. 8. Probability Distribution Functions (PDFs) of gas density (top row) and temperature (bottom row) in the ISM region, weighted by mass for a stack of the simulated galaxies from 300 Myr to 1000 Myr (solid curves) in the ISM regions. The translucent curves are the PDFs of individual snapshots with a spacing of 10 Myr. The x -axes of the T -PDFs are inverted to plot the hot, di ff use gas on the left and the cold, dense gas on the right as the n -PDF. The dotted curves show the stacked results for the low-resolution simulation. The first column shows the PDFs of all the gas in the disk, while the second to fourth columns show the mass distribution of gas for di ff erent galaxies weighted by H ii (purple), H i (orange), and H2 (blue) fractions. The dashed vertical lines on the T -PDFs indicate the boundaries of the CNM, WNM, WIM, and HIM phases, whose corresponding mass and volume fractions in the ISM have been presented in Fig. 6. The ISM in three galaxies with di ff erent metallicities and feedback channels are all dominated by the warm ( ∼ 10 4 K) di ff use ( ∼ 1 cm -3 ) H i gas. Radiative feedback leads to significantly less buildup of H2 in metal-poor environments.  \n<!-- image -->  \nmation of H2 can be e ff ective in the high-temperature di ff use gas with a non-negligible abundance of electrons (Bialy & Sternberg 2019), which should be suppressed in the presence of the ISRF.', '3.4.3. Interstellar radiation field': 'In previous sections, we have seen that the existence of FUV ISRF can significantly change the phase structure of the ISM. In Fig. 9, we present the radial distribution of the FUV ISRF  \nFig. 9. 5.8-13.6 eV ISRF normalized to the Habing (1968) unit as a function of R in the ISM regions for the 0 . 02Z ⊙ / full and 0 . 2Z ⊙ / full galaxies. The solid curves are the stacked results from 300 Myr to 1000 Myr, while the translucent curves are the results of individual snapshots with a spacing of 20 Myr. The dotted curves present the stacked results of the corresponding low-resolution galaxies. The grey dashed line shows the value of the background FUV field of G 0 , UVB = 0 . 00 324.  \n<!-- image -->  \nwithin the ISM regions of the two galaxies with radiative feedback, 0 . 02Z ⊙ / full and 0 . 2Z ⊙ / full. The ISFR energy density is described by the conventional parameter G 0 = u 5.8-13.6 eV / 5 . 29 × 10 -14 erg cm -3 (Habing 1968). The solid curves represent the average ISRF by stacking the snapshots from 300 Myr to 1000 Myr, while the faint curves represent individual snapshots taken with intervals of 20 Myr. In both galaxies, the ISRF decreases gradually from the inner to outer regions, following an approximately exponential decay profile. In the 0 . 2Z ⊙ / full galaxy, the stacked ISRF peaks at around ∼ 0 . 1 G 0 in the central area and decreases exponentially to the background level at R ≃ 2 kpc. The ISRF in 0 . 02Z ⊙ / full is roughly one-tenth of that in the 0 . 2Z ⊙ / full galaxy and decreases from the center to the outer regions in a similar manner (for comparison, the SFR is ∼ 6% of that in 0 . 2Z ⊙ / full). These results show that the ISRF is proportional to the metallicity of the dwarf galaxy. This correlation arises because the FUV radiation originates from young OB stars, making its intensity proportional to the SFR. As reported in section 3.2, the global SFR also shows a positive correlation to the metallicity. However, the increase of ISRF with metalicity is slightly slower than that of SFR, because of the competition between higher SFR and higher dust opacity, along with the reduced FUV and LW photon production rates of metal-rich stars (see Fig. 1). The ISRFs of individual snapshots display an exponential decay trend as well, interspersed with multiple spikes indicating local enhancements in FUV intensity near young massive stars and clusters. At any given radius, there can be temporal variations in the ISRF of more than two or even three orders of magnitude across di ff erent snapshots.  \nThe stacked ISRF in the low-resolution galaxies is also plotted with dotted curves, which show minor di ff erences from the high-resolution results for the 0 . 2Z ⊙ galaxies and the inner part of the 0 . 02Z ⊙ . However, the stacked ISRF in the  \n0  \nFig. 10. Histogram of the hydrogen number density of the host gas cell where the SNe injected their explosion energy. The notations of simulations are identical to those in Fig. 2.  \n<!-- image -->  \n0 . 02Z ⊙ / full / low galaxy appears to be two to four times lower than that of 0 . 02Z ⊙ / full in the outer part ( R ≳ 700 pc, roughly half of the scale length of gas disk). This suggests that the low resolution can suppress star formation in the low gas surface density regions.', '3.4.4. ISM densities around SN explosions': 'One crucial function of early radiative feedback is to disperse the dense star-forming ISM through radiation pressure and the expansion of H ii regions. This process allows SNe to explode in less dense environments, enhancing the e ffi ciency of coupling the explosion energy to the gas. In Fig. 10, we present the distribution of the densities of the gas cells where the SNe released their explosion energy.  \nIn the case of the two full feedback galaxies, most of the SNe exploded in n H < 1 cm -3 di ff use gas because of the preprocessing e ff ect of radiation. The explosion density distributions of these two runs are similar at the high-density end, but for 0 . 2Z ⊙ / full, there are more SNe exploding in n H < 0 . 01, cm -3 low-density gas. This can be understood by the scenario that the radiative feedback in these two simulations rarefies the ISM to a similar degree, so the explosion densities of the first SN in star clusters follow a similar distribution. On the other hand, in 0 . 2Z ⊙ / full, there are more clustered SNe because of its higher SFR, while most star clusters in 0 . 02Z ⊙ / full host only one SN. In a cluster with multiple SNe, the first several SNe will create a hot bubble, allowing the subsequent SNe to explode within this low-density bubble formed by earlier explosions. Consequently, this results in a higher fraction of SNe exploding in low-density environments in 0 . 2Z ⊙ / full.  \nIn the absence of radiative feedback, the explosion density distribution of 0 . 02Z ⊙ / noRT is significantly broader, reaching up to 10 4 cm -3 . Although the main sequence winds can still provide early feedback, their ability to expel cold dense gas is limited by their low luminosities and e ffi cient cooling (Haid et al. 2018; Lancaster et al. 2021a,b). Thus, the first SN of a cluster is likely to occur at a high density and be ine ffi cient in providing feedback and driving outflow (see Section 3.6 for the analysis of outflow properties).  \nFig. 11. Star formation surface density Σ SFR as a function of total (left), thermal (middle), and turbulence (right) pressure of the ISM measured at the galactic mid-plane. The Σ SFR measurement is the same as in Fig. 3. The black solid line in the left panel presents the relation given by Ostriker &Kim(2022) based on simulations of solar metallicity ISM boxes, while the green and red dashed lines present the metallicity-dependent relation given by Kim et al. (2024). We noticed that the numerical experiments performed by Kim et al. (2024) only extend to Z = 0 . 1 Z ⊙ . Therefore, the red dashed line is obtained by extrapolating their fitting function to Z = 0 . 02 Z ⊙ . The full feedback simulations follows this metallicity-dependent superlinear relation, which exhibit steeper slopes and lower intercepts compared to 0.02Z ⊙ / noRT.  \n<!-- image -->', '3.4.5. Interstellar medium pressure modulated by feedback': 'Star formation activities and ISM properties interplay with each other. Stellar feedback injects energy and momentum into the ISM, sustaining its turbulent and thermal pressures, whereas the local deficit of pressure due to cooling leads to the gravitational collapse of the ISM and star formation. This concept is encapsulated in the pressure-regulated, feedback-modulated (PRFM) theory of star-forming ISM (Ostriker et al. 2010; Ostriker & Shetty 2011; Ostriker & Kim 2022).  \nIn vertical dynamical equilibrium of the galactic disk, the weight of the ISM per unit area, W , is balanced by the total ISM pressure, P tot, at the galactic mid-plane. The ISM pressure is sustained by star formation activities that injected momentum and energy through feedback. To satisfy the vertical energy and momentum equations simultaneously, Σ SFR should be proportional to W , and consequently, to P tot. Therefore, the correlation between ISM pressure and star formation surface density serves as a valuable diagnostic for stellar feedback models.  \nIn Fig. 11, we present the Σ SFR as a function of the ISM pressure measured at the galactic mid-plane. Here, the total pressure P tot = P th + P turb, the thermal pressure P th = ρ c 2 s / (5 / 3), and the turbulent pressure P turb = ρσ 2 z are plotted separately. In the two full feedback galaxies, 0 . 02Z ⊙ / full and 0 . 2Z ⊙ / full, the ISM pressure is dominated by the thermal pressure which spans a similar range of 10 2 . 5 Kcm -3 < P / k B < 10 3 . 5 Kcm -3 . The turbulent pressure in these galaxies shows a wider variation, with P turb being overall lower in 0 . 02Z ⊙ / full compared to 0 . 2Z ⊙ / full. The similarity in thermal pressure is due to the prevalence of the warm di ff use ISM, maintaining nearly identical thermal conditions when P / k B < 10 4 Kcm -3 . The lower turbulent pressure in 0 . 02Z ⊙ / full results from the reduced injection of turbulence attributed to its lower SFR. On the other hand, 0 . 02Z ⊙ / noRT exhibits a slightly lower thermal pressure due to the lack of heating by ISRF, while its turbulent pressure is significantly higher compared to 0 . 02Z ⊙ / full, attributed to its enhanced SFR.  \nNow we turn our attention to the correlation between pressure and SFR surface density. The left panel of Fig. 11 shows that the two full feedback simulations roughly follow the superlinear relation between P tot and Σ SFR given by Kim et al. (2024) in their Fig. 10, which exhibit steeper slopes and lower intercepts com-  \nared to 0 . 02Z ⊙ / noRT. The 0 . 02Z ⊙ / noRT simulation roughly follows the relation given by equation (26a) in Ostriker & Kim (2022). We noticed that the numerical experiments performed by Kim et al. (2024) only extend to Z = 0 . 1 Z ⊙ and we extrapolate their metallicity dependence of ∝ Z 0 . 30 to Z = 0 . 02 Z ⊙ to make a comparison with our 0.02Z ⊙ / full galaxy, showing that their fitting functions might still be applicable at metallicity lower than 0 . 1 Z ⊙ . The di ff erence in the intercept between the noRT and full models can be understood qualitatively from two perspectives. Firstly, when considering the same P tot, or in other words, the ISM weight, the full feedback model exhibits a lower SFR due to the impact of radiative feedback in reducing SFE. On the other hand, with the same Σ SFR, the full feedback model will yield more feedback as a result of the presence of radiative feedback and its preprocessing impact on the ISM. The di ff erence in slopes will be explored in a separate study, which will cover a broader dynamical range of P tot and Σ SFR, and will conduct a more detailed quantitative analysis.  \nThe feedback yields ( P / Σ SFR) of the two full feedback galaxies also have di ff erences. The feedback yields for the total pressure are comparable when Σ SFR is high ( ≳ 10 -4 . 5 M ⊙ yr -1 kpc -2 ), but the 0 . 2Z ⊙ / full yields lower pressure in the low Σ SFR regions. This is because the total pressure is dominated by the thermal component. In the case of 0 . 02Z ⊙ / full, feedback is expected to yield higher thermal pressure compared to 0 . 2Z ⊙ / full, especially at lower Σ SFR. This phenomenon can be attributed to ine ffi cient cooling in metal-poor environments, particularly in regions with low Σ SFR where the gas surface density is low, making the gas more susceptible to heating and less prone to cooling. In contrast, feedback in 0 . 02Z ⊙ / full will yield a lower turbulent pressure than in 0 . 2Z ⊙ / full, suggesting a less efficient coupling of feedback to ISM in metal-poor environments. This could be due to the fact that high thermal pressure in lowmetallicity environments results in a less e ffi cient momentum coupling and turbulence driving from stellar feedback (e.g. Haid et al. 2018). The interplay between thermal and turbulent pressure then sets the total ISM pressure which equilibrates the ISM weight and regulates the star formation.  \nFig. 12. Cluster initial mass functions of the simulated galaxies. The solid curves are the results of the high-resolution simulations, while the dotted curves are the results of the corresponding low-resolution galaxies. The grey dashed line presents the observed power-law relation with a slope of -2. Radiative feedback steepens the mass functions by rapidly dispersing the clouds and shutting o ff the further star formation.  \n<!-- image -->', '3.5. The formation and properties of star clusters': "To study the properties of star clusters in our simulations, we perform a four-dimensional 'Friends of Friends' (4D-FoF) analysis on the newly formed stars in the x formt form space, where x form and t form are the formation coordinates and formation time of star particles recorded in the simulation snapshots. Particles with physical connection are linked as a group using a linking length l link and a linking time t link. We choose l link = 5 pc and t link = 5 Myr motivated by the observed size and age spread of star clusters. We found that tuning l link or (and) t link by a factor of two has minor e ff ects on the cluster mass functions and age spreads, while larger or smaller choices result in either overlinking or underlinking. We only keep the identified groups that have at least 35 members as star clusters (Lada & Lada 2003). In the case of low-resolution simulations, we set a minimum membership of 4 to achieve a similar minimum cluster mass to that of high-resolution simulations.", '3.5.1. Cluster initial mass functions': 'In Fig. 12, we present the cluster initial mass functions (CIMFs) of our simulated galaxies. The initial mass of each cluster is defined as the sum of the initial mass of its member star particles. In the 0 . 2Z ⊙ / full galaxies, the CIMF follows a power-law distribution of d N / d M cl ∝ M -α cl with a slope α = 2 . 07, which is in agreement with the observed -2 slope (e.g. Portegies Zwart et al. 2010; Krumholz et al. 2019). The CIMF of the 0 . 02Z ⊙ / full galaxies has a slightly steeper slope of -2 . 40, although the uncertainty of the slope is large due to low number statistics. Turning o ff the radiative feedback results in a significantly shallower power-law slope of -1.63, indicating that more massive clusters are formed. Thus, we conclude that radiative feedback reduces the clustering of star formation by decreasing the number of  \n]  \nr  \ny  \nM  \n[  \nt  \nFig. 13. Age spreads of the clusters ( t 90 -10) as the function of the cluster initial mass. The data with errorbars present the median age spreads of the [0,20), [20,40), [40,60), [60,80), [80,100] cluster initial mass percentile bins, and the vertical errorbars show the 16 and 84 percentiles of the age spreads in each bin. More massive clusters exhibit larger age spreads, suggesting that the final mass of the cluster is influenced by the duration of star formation within the cloud when there is no disruptive stellar feedback.  \n<!-- image -->  \nmassive clusters and shaping the CIMF steeper. This indicates that radiative feedback disperses clouds quickly and e ffi ciently.  \nThe power-law slopes for the low-resolution simulations are -2 . 16 for 0 . 2Z ⊙ / full / low, -2 . 34 for 0 . 02Z ⊙ / full / low, and -1 . 80 for 0 . 02Z ⊙ / noRT / low, respectively. Despite a minor deviation from the high-resolution values, our qualitative findings regarding the shapes of their mass distributions remain robust to numerical resolution.', '3.5.2. Age spreads and feedback timing': "Now, we study the cluster formation time based on the age spreads of the member stars. We use the age spreads as a proxy to study the time-scale of feedback dispersing the star-forming clouds. We quantify the age spread of a cluster by estimating the duration needed for a specific percentage of member stars to be formed within the cluster. Since the mass of each star particle (inherited from its progenitor gas cell, not M fb) is almost the same, it describes a mass-weighted cluster formation timescale. In Fig. 13, we plot the age spreads of the cluster when 10 per cent to 90 per cent of its final members are formed ( t 90 -10) as a function of the cluster initial mass. As can be seen, the age spread shows a clear positive correlation with the cluster initial mass, where the more massive clusters exhibit larger age spreads, especially for the 0 . 02Z ⊙ / noRT. This correlation suggests that the final mass of the cluster is influenced by the duration of star formation within the cloud in the absence of disruptive stellar feedback.  \nTo quantitatively analyze the correlation between feedback and cluster growth, we study the formation time of the first massive star in the cluster ( t fms). In Fig. 14, we present the duration between the formation of the first and last member star of a clus-  \nFig. 14. Time between the formation of the first and last member star of a cluster ( t 100 -0) as the function of the time of the first massive star in the cluster ( t fms). The black dashed curves present the line of t 100 -0 = t fms + { 0 Myr , 1 Myr , τ⋆ (100 M ⊙ ) , τ⋆ (20 M ⊙ ) } , where τ⋆ (100 M ⊙ ) and τ⋆ (20 M ⊙ ) are the lifetimes for 100 M ⊙ and 20 M ⊙ massive stars. The time between t 100 and t fms are smaller than 1 Myr for most clusters in the full feedback simulations, much shorter than the time required for the first SN explosion.  \n<!-- image -->  \nter ( t 100 -0) as a function of t fms. The timing at which SN feedback occurs is plotted with black dashed curves. Among these curves, the ' + τ⋆ (100M ⊙ )' curve corresponds to the minimum explosion time of the first potential SN in the cluster, given that 100 M ⊙ represents the maximum stellar mass in our model. As shown by the results of 0 . 02Z ⊙ / full and 0 . 2Z ⊙ / full, the time differences between t 100 and t fms are smaller than 1 Myr for most clusters in the full feedback simulations, much shorter than the time required for the first SN explosion. The clusters cease to grow rapidly after the formation of the first massive stars, suggesting that the clouds are dispersed rapidly by the pre-SN radiative feedback. Without radiative feedback, the clouds can remain in their star-forming phase for longer periods. While the stellar winds can also contribute to early feedback, its dynamic influence is relatively ine ff ective. Therefore, the clusters grow over long periods and accumulate higher final masses.  \nIn summary, in low-surface density, low-metallicity environments such as our dwarf galaxies, the mass growth period and the final mass of low-mass ( < 10 4 M ⊙ ) young star clusters are closely correlated with the timing of feedback emergence. With e ff ective radiative feedback, the appearance of the first massive star nearly coincides with the dispersal of its birth molecular cloud, leading to the inhibition of further star formation. SN feedback typically does not contribute to cloud dispersal. The absence of radiative feedback allows the cloud to exist longer and form new stars, and SNe also play a role in disrupting clouds in this case.", '3.6. Outflow properties': 'In previous subsections, we show that radiative feedback and metallicity largely a ff ect star formation activities, the ISM struc-  \nture, and the properties of star clusters. Now we investigate how they shape the large-scale outflow properties in the disk-halo interface of the simulated dwarf galaxies. We compute the outflow rates of a physical quantity q (i.e., gas mass ˙ M out, metal mass ˙ M Z,out, momentum ˙ p out, energy ˙ E out) as the sum of qk over all gas particles within a slab of thickness ∆ L centered at height | z |  \n˙ q out = 1 ∆ L X k v k · n qk , ( v k · n > 0) , (17)  \nwhere v k is the velocity of the gas cell k and n is the outward norm of the mid-plane. We measure the outflow rates at | z | = 1 kpc to capture the outflow properties at the launching interface, and we set ∆ L = 100 pc.  \nWe use loading factors to describe how star formation activities drive the galactic outflow. We define the mass-loading ( η M ), metal-loading ( η Z ), momentum-loading ( η p ), and energyloading ( η E ) factors as  \nη M = ˙ M out ˙ M ⋆ , η Z = ˙ M Z,out Z ISM ˙ M ⋆ , η p = ˙ p out p ej ˙ N SN , η E = ˙ E out E SN ˙ N SN , (18)  \nrespectively, where the SFR ˙ M ⋆ is time-averaged over 10 Myr, Z ISM is the metallicity of the ISM region, ˙ N SN = ˙ M ⋆/ (100 M ⊙ ) is the SN rate assuming a canonical IMF, p ej = 3 × 10 4 M ⊙ kms -1 and E SN = 10 51 erg are the momentum carried by ejecta and injected energy of a typical SN.  \nIn the left column of Fig. 15, we present the outflow rates of gas mass, metal mass, momentum, and energy in 10 Myr steps for our simulated galaxies. The low-resolution results are presented with dotted lines, which are similar to the high-resolution results. The 0 . 02Z ⊙ / noRT run has the highest outflow rates for all quantities due to its highest SFR. The outflow rates of the two full feedback galaxies show evident temporal fluctuations, especially the 0 . 02Z ⊙ / full ones, due to their bursty star formation. On the contrary, the 0 . 02Z ⊙ / noRT results present the smallest fluctuations due to its more continuous star formation.  \nIn the right column of Fig. 15, we present the loading factors defined by equation (18). The two full feedback galaxies exhibit comparable (less than 1 dex di ff erence) loading factors for mass, metal, and energy, despite di ff erences in their metallicities and SFRs. The momentum-loading factor of 0 . 02Z ⊙ / full is slightly higher than that of 0 . 2Z ⊙ / full. The similarity in these loading factors is attributed to the stellar feedback mechanisms driving the outflow. The early radiative feedback preprocesses the ISM and reduces the clustering of stars so that sparse SNe explode in less dense gas, though the details of feedback budgets and cooling rates are di ff erent in these two galaxies. The mass- and metalloading factors of 0 . 02Z ⊙ / full and 0 . 2Z ⊙ / full fluctuate rapidly between 10 and 500, with a median of η M ≈ 46, η Z ≈ 45 for 0 . 02Z ⊙ / full and η M ≈ 22 and η Z ≈ 25 for 0 . 2Z ⊙ / full, respectively. Such high loading factors indicate that the gas and metals expelled from the ISM far exceed the mass consumed in star formation in dwarf galaxies. The energy-loading factors of these two runs are similar as well and vary between 0 . 01 and 1 with a median of ∼ 0 . 1, which are significantly smaller than the momentum-loading factors that vary mostly between 1 and 100. This di ff erence can be explained by the fact that the momentum of the SN ejecta is boosted by converting the SN thermal energy to the kinetic energy of the gas through P d V work, while the total energy dissipates during the outflow process.  \nRemarkably, the mass-loading factor of 0 . 02Z ⊙ / noRT is lower than that of the full feedback simulations by almost an  \nFig. 15. Outflow properties of the simulated dwarf galaxies. Left: the outflow rate of gas mass, metal mass, momentum, and energy measured at | z | = 1 kpc. Right: the mass-loading, metal-loading, momentum-loading, and energy-loading factors computed as outflow rate per locking rate as defined in equation (18). Dashed horizontal lines indicate the median values of outflow rates and loading factors in the high-resolution simulations after 100 Myr (outside grey shaded region).  \n<!-- image -->  \n<!-- image -->  \nE  \norder of magnitude with a median of 7, while the metal- and energy- loading factors are similar to those of the other two runs. The reduction of mass and momentum loading factors can be intuitively explained by the scenario in which more SNe explode in dense gas without early radiative feedback, making it more challenging for them to drive outflows (e.g., H17). However, this finding contradicts the recent results of Smith et al. (2021), who suggest that photoionization feedback could suppress galactic outflows by reducing the clustering of SNe. Indeed, radiative feedback prohibits the formation of massive clusters so that de-clustering the SNe in our simulations as well (Section 3.5). However, this e ff ect is mitigated in our simulations due to the low gas-surface density (3 M ⊙ pc -2 at the center) of our galaxies, which limits the ability of clustered SNe to drive outflows. The mass-loading factor becomes saturated once the gas has been cleared by multiple SNe, causing subsequent SNe in the same cluster to directly eject their material into the CGM. Moreover, all our galaxies disfavour the formation of large clumps and mas-  \nve clusters due to the relatively low gas surface density and gas fraction of our initial condition (see Renaud et al. 2021, for the impact of gas fraction). Even in 0 . 02Z ⊙ / noRT, only 6 per cent of the clusters are more massive than 1000 M ⊙ which allows them to host ten SNe. In large clusters, clustered SNe can occur within the low-density bubble created by a previous SN, leading to more e ff ective feedback despite the lack of radiative feedback. In contrast, smaller clusters with a few SNe may struggle to provide substantial feedback. Thus, the ability to form su ffi -cient large clusters also determines the e ff ectiveness of SN explosion in driving galactic outflows in the noRT case.  \nThe loading factors in the low-resolution simulations show a good convergence with the high-resolution ones. Variations in the mass-, metal-, and momentum-loading factors remain below a factor of 1.5, while the energy-loading factors are systematically lower than the high-resolution ones by a factor of up to 2.  \n4', '4. Discussion': '4.1. Effects of metallicity and radiative feedback on simulated dwarf galaxies  \nIn this work, we simulated two galaxies alongside our fiducial simulation, one with ten times higher metallicity (0 . 2Z ⊙ / full) and another with radiative feedback turned o ff (0 . 02Z ⊙ / noRT), to study the e ff ects of metallicity and radiative feedback.  \nThe ten-times-higher metallicity results in a roughly tentimes-higher metal cooling rate in T > 10 5 K hot gas and T < 10 4 Kwarm and cold gas. The gas heated by stellar feedback can rapidly cool down, leading to a new cycle of star formation (Semenov et al. 2017). As a result, the SFR increases by a factor of 17 due to the availability of more cold gas fuel for star formation, and the intensity of ISRF also increases correspondingly by a factor of 10. The increase of ISRF is slightly lower than that of SFR. This can be explained by the decrease of FUV and LW photon production rates with metallicity and the increase of dust absorption with metallicity. Nevertheless, both the KS and extended KS relations of 0 . 02Z ⊙ / full and 0 . 2Z ⊙ / full are consistent with observed dwarf galaxies (e.g. Roychowdhury et al. 2015) and exhibit superlinear relationships with similar slopes. The results of 0 . 2Z ⊙ / full roughly follow the empirical relations given by Elmegreen (2015) and Shi et al. (2018), while those of 0 . 02Z ⊙ / full show a smaller intercept compared with these relations, suggesting a longer dynamical time-scale for star formation in extremely metal-poor environments. The metal-poor and metal-rich galaxies share several properties in common, including the rapid dispersion of clouds, a power-law CIMF with a steep slope, and the loading factors. This similarity arises from the fact that radiative feedback operates similarly on a cluster scale, and the di ff erence in ionizing photon productions between 0 . 02 Z ⊙ and 0 . 2 Z ⊙ is within a factor of three. On the other hand, more massive clusters with larger age spreads can be formed in 0 . 2Z ⊙ / full because this galaxy can form larger molecular clouds that require more energetic feedback to disperse.  \nA key feature of radiative feedback is that it acts immediately following the birth of massive stars, without any delay. On the cloud scale, radiative feedback can rapidly disrupt the star-forming clouds through photoionization and direct radiation pressure, thereby halting further star formation within approximately 1 Myr after the emergence of massive stars (see Section 3.5.2). Thus, a small portion of the star-forming gas is converted into stars before being dispersed by radiative feedback. The growth of cluster masses is also prevented after radiative feedback is initiated, leading to a steep slope of ≲ -2 in the CIMF. However, this situation may not apply in denser highredshift environments, where FUV radiation from nearby massive stars can even enhance the formation of massive clusters by delaying further star formation and enabling the clouds to accumulate more gas (e.g. Garcia et al. 2023; Sugimura et al. 2024). Future investigations will aim to understand the impact of radiative feedback on star cluster formation in various environments.  \nOn the galaxy scale, the local and long-range photoheating processes remove the cold fuel for star formation and prevent the cooling of warm ISM. These combined local and global e ff ects of radiative feedback consequently reduce the SFR in our dwarf galaxy by a factor of 92 when compared to the noRT model. The inclusion of radiative feedback naturally regulates the star formation process in consistency with empirical relationships. In terms of the galactic outflows, the outflow rates are suppressed in models with radiative feedback because of the lower SFR. On the other hand, radiative feedback preprocesses the environments of SN explosion, causing most SNe explode in n H < 1 cm -3 low-  \n0  \nFig. 16. Stellar metallicity distribution of the simulation without any early feedback mechanism (black curve, 0 . 02Z ⊙ / noESF) in comparison with 0 . 02Z ⊙ / noRT and 0 . 02Z ⊙ / full. Stars with extremely high metallicity (log ( Z / Z ⊙ ) > 1) can be formed from the unmixed ejecta in 0 . 02Z ⊙ / noESF due to the absence of early feedback.  \n<!-- image -->  \ndensity gas. This e ff ect increases the occurrence of e ffi cient SN explosions, which in turn increases the mass loading factor of our simulated dwarf galaxy.', '4.2. Necessity of early feedback': "We emphasize that early feedback is indispensable for highresolution multiphase ISM simulations. The maximum gas density of the simulation is constrained by the Jeans criterion ( m gas > f 3 J,s M J, see Section 2.3), indicating that with higher mass resolutions, stars can be formed in significantly denser gas, which scales with m 2 gas . In the absence of e ff ective early feedback, the explosion density of SNe is linked to the density at which they are formed. For example, if a SN explodes in gas with a density of 10 6 cm -3 , the cool radius and cooling time of its ST phase are only 0 . 07 pc and 22 yr (equation (7) and (8) from Kim & Ostriker 2015). The hot SN remnant can only occupy a very small volume for a very short time before rapidly cooling down to cold gas. Although the SN can still inject a substantial amount of momentum ( p final ∝ n -1 / 7 ) as long as the ST phase is resolved, it does not heat the gas to the hot and warm phases. Consequently, the SN ejecta may be misidentified as star-forming gas because of its cold and dense nature. Stars can therefore be sampled from the unmixed ejecta and result in the formation of star particles with extremely high metallicities.  \nTo demonstrate this issue, we run an additional simulation without any early feedback mechanism (i.e., SN only), and present the distribution of stellar metallicity in Fig. 16. The SNonly simulation is denoted as 0 . 02Z ⊙ / noESF, in comparison to 0 . 02Z ⊙ / noRT and 0 . 02Z ⊙ / full. Stars with extremely high metallicity (log ( Z / Z ⊙ ) > 1) can be formed from the unmixed ejecta in 0 . 02Z ⊙ / noESF due to the absence of early feedback. Such stars have never been observed so far. In contrast, the stellar metallicity distribution in 0 . 02Z ⊙ / full exhibits a spread of ∼ 0 . 2 dex. While 0 . 02Z ⊙ / noRT shows a broad range of stellar metallicities extending to the solar value, the inclusion of wind feedback prevents the formation of excessively metal-rich stars. It is important to note that these phenomena are not a result of numerical over-cooling issues but rather the 'physical' results of SN exploded in extremely dense gas. Therefore, we must turn at least the wind feedback on to avoid SNe occurring in very dense gas.", '4.3. Comparison to previous work': "The initial conditions of our simulations are intentionally set to be similar to that in H17 and L23 for a direct comparison. The SFRs of 0 . 2Z ⊙ / full and 0 . 02Z ⊙ / full show good agreement with the PE-PI-SN galaxy in H17 and the 'extended' galaxy in L23. This result is expected since all these simulations incorporate detailed stellar feedback, particularly the pre-SN early feedback, and accurately model the thermodynamics of the multiphase ISM. Consistent with simulations with explicit RT (e.g., Emerick et al. 2018; E19; Andersson et al. 2024) or approximate radiative feedback (e.g., H17; Smith et al. 2021), our simulations also demonstrate that the inclusion of radiative feedback alongside SN feedback leads to a reduction in the SFR. However, in Emerick et al. (2018), the radiative feedback only decreased the SFR by about a factor of 5 compared to their 'No RT' run, whereas in our simulations, this reduction is 92. We note that Emerick et al. (2018) simulated a much smaller system than ours. In their 'No RT' run, the SFR experiences a significant decline after the initial 50 Myr and does not reach a selfregulating state. It is unknown if the SFR could further increase when more gas is available, leading to better agreement with our results.  \nThe warm H i gas dominates the ISM in all three simulated dwarf galaxies, which is consistent with previous simulations of dwarf galaxies (e.g., Hu et al. 2016; H17; E19). Hot ionized gas also fills a non-negligible ISM volume (up to ∼ 0 . 1) in our 0 . 2Z ⊙ / full galaxies, as reported by H17. In contrast, the volume fraction of HIM is only ∼ 10 -3 in 0 . 02Z ⊙ / full due to its low SFR. In line with the findings of Hu et al. (2021), we found that there is a substantial fraction of H2 mass in warm di ff use gas with a low H2 fraction in low-metallicity environments.  \nThe ISRF of our 0 . 2Z ⊙ / full presents a similar radial profile as that of PE-PI-SN in H17. The intensity of our ISRF is slightly lower than that of H17 at a given radius. This can be attributed to the slightly higher metallicity of our galaxy, which results in more dust absorption. The significant temporal fluctuations in our ISRF also align with the findings of E19, who similarly reported a variation exceeding two orders of magnitude.  \nSimilar to Andersson et al. (2024), we found that the masses and time spreads of star clusters are sensitive to the timing of the onset of radiative feedback. Andersson et al. (2024) also found that radiative feedback results in a steeper power-law slope of the CIMF. However, they report a slope around -3 for their 'SNe + winds + radiation' model. Despite this discrepancy, our work again underscores the importance of incorporating the early feedback in galaxy formation models from the perspective of star cluster formation.  \nThe mass-loading factor of our 0 . 02Z ⊙ / full galaxy varies in the order of 10 -500, which is similar to the results of the 'extended' galaxy in L23. Emerick et al. (2018) also reported a comparable range for the mass-loading factor in their fiducial model. We get enhanced mass-loading factors with radiative feedback turned on, and we attribute this result to the early radiative feedback that enable SNe to explode in less dense gas with n H < 1 cm -3 so that they are more e ffi cient to drive galactic outflow. This pre-pocessing e ff ect of radiative feedback on the SN explosion site has been noticed in simulations of various environments (e.g., H17, Kannan et al. 2020a, Smith et al. 2021, Rathjen et al. 2021). Notably, the 'shortrad' model in Emerick et al. (2018), which modelled short-range radiative feedback but lacked the long-range e ff ects, exhibited a lower mass-loading factor ranging from 1 to 10. This result emphasizes the importance of explicit RT as we have integrated into the RIGEL model.  \nAs we mentioned in Section 3.6, our enhanced mass-loading factor in the full feedback simulations is opposite to the finding by Smith et al. (2021). We explain this di ff erence with the low gas surface density and (relatively) low gas fraction of our initial condition, leading to star clusters that have only a few SNe, even in the noRT simulation. Simulations of higher gas surface density environments might show di ff erent behavior. Consistent with our results, the simulations of a smaller galaxy in Emerick et al. (2018) also reported a significantly reduced mass loading factor for their 'No RT' run. However, they halted this run after ∼ 100 Myr evolution so the conclusion may still be unclear. Moreover, Smith et al. (2021) only models the short-range effects of photoionization using the Strömgren-type approximation, which can reduce the loading factor as reported by Emerick et al. (2018). We will discuss the dependence of the outflow properties on the gas fraction and clustering properties in an accompanying paper.", '4.4. Limitations and future perspectives': "We note that there are a few caveats and limitations of this work. Firstly, since the resolution of our simulations is insu ffi cient to resolve the formation and ongoing accretion of individual protostars, we use a sub-grid model outlined in Section 2.3 to model star formation following the Schmidt (1959) law. Although our model successfully reproduces the observed star formation relations in dwarf galaxies as shown in Section 3.2, it does not provide insights into the star formation properties such as the spatial and temporal preferences of star formation across varying stellar mass ranges and the fundamental physics behind the IMF. As a result, though the global properties like SFR are insensitive (Hopkins et al. 2011), the properties within star clusters, such as the environmental densities of star formation, compactness and boundness of star clusters, and cluster formation e ffi ciency, can vary significantly with the underlying sub-grid model (Hislop et al. 2022). Moreover, the IMF can be more topheavy than the canonical one in metal-poor dwarf galaxies (e.g. Marks et al. 2012; Gennaro et al. 2018), thereby increasing the e ffi ciency of stellar feedback in regulating star formation (e.g. Gutcke & Springel 2019; Prgomet et al. 2022). Although the incorporation of a variable IMF is relatively straightforward, the relationship between the slope of the IMF and physical properties such as metallicity remains unclear. On the other hand, a variety of observational clues also suggest a relation between the highest initial stellar mass and the cluster mass (Larson 2003; Yan et al. 2023). As a result, realistic low-mass clusters may not be able to host very massive stars stochastically sampled from the IMF. More sophisticated star formation models can help us capture this observed relation at the cost of computational e ffi -ciency (e.g., Hirai et al. 2021; L23).  \nRecent observations suggest that the classical stellar wind model we adopted can overestimate wind luminosities (Smith 2014). Although there are still large theoretical and observational uncertainties, we have compiled a table of fitting parameters based on the weak-wind model following Grudi'c et al. (2021) for future study. We also do not model post-mainsequence winds because i) in our metal-poor dwarfs, only a handful number (order of 10) of ≳ 20 M ⊙ stars can be formed during 1 Gyr, which have a Wolf-Rayet (WR) phase; ii) for these stars, though their mass-loss rates are enhanced by a factor of ∼ 10, they have similar wind luminosities to the main-sequence winds (e.g., Meynet & Maeder 2005); iii) the duration of the WR phase is short ( ≲ 1 Myr, LC18). Thus, we argue that the post-main-sequence winds have negligible e ff ects on the global  \nproperties of our simulated galaxies. We assume that massive stars release unenriched gas by main-sequence winds and release yields by post-main-sequence winds together with SN ejecta. This assumption is valid for the non-rotating massive stars in the LC18 model. However, fast-rotating stars with e ffi cient mixing and strong winds can provide continuous yield input throughout their lifetime. Moreover, we do not include binary stars in our model, while feedback from binaries can be very di ff erent from that of single stars in terms of production rates of ionizing photons and metal yields (e.g. Götberg et al. 2017; Secunda et al. 2020; Tsai et al. 2023; Yates et al. 2024). The enrichment contributions of fast-rotating stars and binary stars are possibly important ingredients in producing multiple populations in clusters (Bastian & Lardo 2018). One ongoing investigation is implementing the post-main-sequence evolution and yields from fast-rotating stars and binaries to study the multiple populations problem.  \nWe assume that all 8 -100 M ⊙ massive stars have a constant SN energy of 10 51 erg. However, the 'explodability' of massive stars and the type of their remnants can be a complicated function of their masses and metallicities (Heger et al. 2003), and the explosion energy can be a spectrum over a wide range (e.g. Ertl et al. 2016; Sukhbold et al. 2016). Gutcke et al. (2021) modelled the SN explosion with variable energy as a function of mass based on Sukhbold et al. (2016). They found that their model yielded a factor of 2 -3 lower mass loading factor compared to the fixed SN energy model. The implementation is straightforward for us, but the metallicity dependence of SN energy spectrum across the entire mass range is still unclear.  \nThe metals di ff use passively following mesh advection, and we do not employ any sub-grid di ff usion model for the unresolved eddies. In this case, the numerical di ff usivity is roughly on the scale of σ t ∆ x , where σ t is the velocity dispersion and ∆ x is the spatial resolution. Recent work by Steinwandel et al. (2024) has shown that the inclusion of a metal di ff usion model has a minor impact on ISM structures and loading factors, while it plays a significant role in the phase structure of the galactic outflows. Investigating metal di ff usion in detail will be a significant focus of the RIGEL project.  \nDust formation and destruction in the gas are not explicitly accounted for in our current model. We assumed a constant dustto-metal ratio of 0.5 based on the MW observations. However, Rémy-Ruyer et al. (2014) suggest a lower dust-to-metal ratio at low metallicity. To improve the accuracy of our dust physics modeling, future iterations of the model will need to explicitly consider processes such as dust formation and destruction, as well as the evolution of dust size and its interaction with radiation (e.g McKinnon et al. 2018, 2021). These advanced features will be integrated into our model in the future.  \nLastly, in this work, we focused on studying an idealized dwarf galaxy at z = 0. Therefore, the conclusions are mostly limited to this type of dwarf galaxy with relatively quiescent star formation. Future works on simulating dwarf galaxies in di ff erent environments, such as high-z gas-rich analogs, merging systems, and ultimately the cosmological zoom-in boxes are needed to investigate the stellar feedback and multiphase ISM in dwarf galaxies across cosmic time. Specifically, we assume a uniform z = 0 UVB in this work, which is appropriate for our non-cosmological isolated dwarf galaxy. However, the spatial and temporal variations of UVB have a significant impact on the low-mass galaxy populations, a ff ecting their abundance, star formation history, physical extent, and metal enrichment (Kim et al. 2023b; Borrow et al. 2023). One immediate application of RIGEL model is to explore the cosmological evolution of the  \nlow-mass reionization survivors ( M halo ≳ 10 9 M ⊙ , Jeon et al. 2017, 2021; Gutcke et al. 2022; Lee et al. 2024) in synergy with the inhomogeneous UVB predicted by the THESAN simulations (Kannan et al. 2022; Smith et al. 2022; Garaldi et al. 2022).", '5. Summary': 'This paper introduces the RIGEL model, a framework for studying self-regulating star formation, multiphase ISM, galactic outflows, and cluster formation in dwarf galaxies and other galactic environments. The primary advancement of RIGEL, in comparison to similar numerical models, lies in its incorporation of explicit RT and a detailed feedback model for massive stars across a broad range of metallicities from zero to solar. We incorporate metallicity-dependent, multi-channel stellar feedback on a starby-star basis, whereby we explicitly sample individual massive stars from the IMF to account for feedback mechanisms such as radiation, stellar winds, and core-collapse SNe. This treatment of feedback, combined with the e ffi cient M1 RT solver and improved heating-cooling model that consider nonequilibrium primordial chemistry and equilibrium abundances of C / O species, enables us to comprehensively model self-regulating star formation processes and the thermodynamics of all ISM phases, spanning from cold molecular gas to hot ionized gas. In this work, we validate our model by simulating the evolution of an isolated dwarf galaxy at gas mass resolutions of 1 M ⊙ and 10 M ⊙ . For our simulated dwarf galaxies, we find that:  \n- 1. The RIGEL model self-consistently produces a multiphase ISM (Fig. 4 and 5) in which the star formation is regulated to a rate consistent with the observed dwarf galaxies (see Section 3.2).\n- 2. The SFRs of our simulated dwarf galaxies show a strong positive correlation with the metallicity. Radiative feedback can reduce the SFR by almost two orders of magnitude (see Fig. 2 and 3) by removing the cold, dense gas fuel available for star formation (Fig. 6). Turning o ff radiative feedback boosts the SFR to an unreasonably high value by about two orders of magnitude.\n- 3. The ISMs in three galaxies with di ff erent metallicities and feedback channels are all dominated by warm (100 K < T < 10 5 K) H i gas in both mass and volume. Specifically, H i gas dominates at densities from n H ∼ 10 -2 . 5 cm -3 to at least 10 3 cm -3 , covering a broad range of temperatures of T ∼ 10 -10 4 K. H ii gas dominates the hot and photoionized gas. H2 is always a minority phase in the 0 . 02Z ⊙ / full galaxy, even at a high density of n H > 10 3 cm -3 (see Fig. 8). With a higher metallicity, 0 . 2Z ⊙ / full is still dominated by H i below 10 3 cm -3 . Turning o ff radiative feedback leads to di ff use H2 formed in warm-hot gas with temperatures and densities extended to 10 6 K and 10 -4 cm -3 in the absence of the LW radiation.\n- 4. The ISRFs of our galaxies have significant temporal variations of two orders of magnitude, but on average they have an exponentially declining radial profile that peaks at the galactic center (see Fig. 9). The average ISFR of 0 . 2Z ⊙ / full is roughly ten times higher than 0 . 02Z ⊙ / full because of its higher SFR.\n- 5. Radiative feedback plays a crucial role in removing the dense gas from the star-forming regions to create rarefied gas cavities where SNe explode (see Fig. 10). Although stellar winds may not be very e ffi cient in dispersing dense gas, they still play an important role in preventing SN explosions in highly dense environments. In Section 4.2, we show that turning o ff  \nall the early feedback results in a rapid cooling of (unmixed) SN ejecta within dense gas, leading to the formation of excessively metal-rich stars.  \n- 6. The ISM pressure of our three galaxies is dominated by the thermal component. The total pressure and star formation surface density follow a superlinear relation. The full feedback galaxies have a lower SFR at a given pressure compared to the noRT one, due to the impact of radiative feedback in reducing SFE. On the other hand, with the same star formation surface density, the full feedback galaxies yield more pressure because of the presence of radiative feedback and its preprocessing e ff ect on the ISM (see Section 3.4.5).\n- 7. The emergence of radiative feedback can rapidly disperse the molecular clouds and suppress further star formation. The age spread of low-mass star clusters is therefore tightly related to the timing of the onset of radiative feedback (Fig. 13). The e ffi cient early feedback prohibits the formation of massive star clusters, shapes the CIMF slope at the high-mass end (Fig. 12), and shrinks the age spread of star clusters to less than 2 Myr (Fig. 14).\n- 8. The mass-loading factors of the galactic outflows in our full feedback galaxies vary between 10 and 500 (Fig. 15). Turning o ff the radiative feedback leads to reduced mass- and momentum-loading factors, while having minor e ff ects on the metal- and energy-loading factors. The relationship between outflow properties and the properties of galactic disks will be explored in an accompanying paper.\n- 9. Finally, the 10 M ⊙ low-resolution simulations present a satisfactory convergence with the 1 M ⊙ fiducial simulations concerning various aspects such as the SFRs, the ISM structures, the ISRF distribution, the cluster mass functions, and the outflow properties. This convergence enables us to investigate the evolution of larger systems over cosmic time with a resolution of 10 M ⊙ with high fidelity.  \nAcknowledgements. Wethank Volker Springel for giving us access to arepo . YD is grateful to Yang Ni and Zhiqiang Yan for useful discussions. HL is supported by the National Key R&D Program of China No. 2023YFB3002502, the National Natural Science Foundation of China under No. 12373006, and the China Manned Space Program through its Space Application System. BL is supported by the Royal Society University Research Fellowship. RK acknowledges support of the Natural Sciences and Engineering Research Council of Canada (NSERC) through a Discovery Grant and a Discovery Launch Supplement, funding reference numbers RGPIN-2024-06222 and DGECR-2024-00144. The simulations of this work were run on the Stampede2 HPC resource at the Texas Advanced Computing Center and the Anvil cluster at Purdue University as part of ACCESS through TG-MCA06N030 and TG-PHY220084. We use python packages N um P y (Harris et al. 2020), S ci P y (Virtanen et al. 2020), astropy (Astropy Collaboration et al. 2013, 2018), and matplotlib (Hunter 2007) to analyze and visualize the simulation data.', 'References': "- Agertz O., Kravtsov A. V., 2015, ApJ, 804, 18\n- Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25 Agertz O., et al., 2020, MNRAS, 491, 1656\n- Andersson E. P., Agertz O., Renaud F., Teyssier R., 2023, MNRAS, 521, 2196\n- Andersson E. P., Mac Low M.-M., Agertz O., Renaud F., Li H., 2024, A&A, 681, A28\n- Applebaum E., Brooks A. M., Quinn T. R., Christensen C. R., 2020, MNRAS, 492, 8  \nAsplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481  \n- Astropy Collaboration et al., 2013, A&A, 558, A33\n- Astropy Collaboration et al., 2018, AJ, 156, 123  \nAubert D., Teyssier R., 2008, MNRAS, 387, 295  \n- Bagla J. S., 2002, Journal of Astrophysics and Astronomy, 23, 185  \nBarnes J., Hut P., 1986, Nature, 324, 446  \nBastian N., Lardo C., 2018, ARA&A, 56, 83  \nBialy S., Sternberg A., 2019, ApJ, 881, 160  \n- Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2846\n- Bigiel F., Leroy A., Walter F., Blitz L., Brinks E., de Blok W. J. G., Madore B.,\n- 2010, AJ, 140, 1194\n- Bisbas T. G., et al., 2015, MNRAS, 453, 1324\n- Bode P., Ostriker J. P., Xu G., 2000, ApJS, 128, 561\n- Borrow J., Kannan R., Garaldi E., Smith A., Vogelsberger M., Pakmor R., Springel V., Hernquist L., 2023, MNRAS, 525, 5932\n- Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni\n- A., 2012, MNRAS, 427, 127\n- Bruzual G., Charlot S., 2003, MNRAS, 344, 1000\n- Burke J. R., Hollenbach D. J., 1983, ApJ, 265, 223\n- Calura F., et al., 2022, MNRAS, 516, 5914\n- Cassisi S., Castellani V., 1993, ApJS, 88, 509\n- Cen R., 1992, ApJS, 78, 341\n- Chabrier G., 2003, PASP, 115, 763\n- Chan T. K., Theuns T., Bower R., Frenk C., 2021, MNRAS, 505, 5784\n- Chevance M., et al., 2022, MNRAS, 509, 272\n- Chevance M., Krumholz M. R., McLeod A. F., Ostriker E. C., Rosolowsky E. W., Sternberg A., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, Astronomical Society of the Pacific Conference Series Vol. 534, Protostars and Planets VII. p. 1 ( arXiv:2203.09570 ), doi:10.48550 / arXiv.2203.09570\n- Conroy C., Wechsler R. H., 2009, ApJ, 696, 620  \nDekel A., Silk J., 1986, ApJ, 303, 39  \n- Deng Y., Li H., Kannan R., Smith A., Vogelsberger M., Bryan G. L., 2024, MNRAS, 527, 478\n- Dobbs C. L., Burkert A., Pringle J. E., 2011, MNRAS, 417, 1318\n- Doherty C. L., Gil-Pons P., Lau H. H. B., Lattanzio J. C., Siess L., Campbell\n- S. W., 2014, MNRAS, 441, 582\n- Draine B. T., 1978, ApJS, 36, 595\n- Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium\n- Eftekhari F. S., et al., 2022, MNRAS, 517, 4714\n- Elmegreen B. G., 2015, ApJ, 814, L30\n- Elmegreen B. G., 2018, ApJ, 854, 16\n- Elmegreen B. G., Hunter D. A., 2015, ApJ, 805, 145\n- Emerick A., Bryan G. L., Mac Low M.-M., 2018, ApJ, 865, L22\n- Emerick A., Bryan G. L., Mac Low M.-M., 2019, MNRAS, 482, 1304\n- Ertl T., Janka H. T., Woosley S. E., Sukhbold T., Ugliano M., 2016, ApJ, 818, 124\n- Faucher-Giguère C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ, 703, 1416\n- Faucher-Giguère C.-A., Quataert E., Hopkins P. F., 2013, MNRAS, 433, 1970\n- Fishlock C. K., Karakas A. I., Lugaro M., Yong D., 2014, ApJ, 797, 44\n- Fujimoto Y., Chevance M., Haydon D. T., Krumholz M. R., Kruijssen J. M. D., 2019, MNRAS, 487, 1717\n- Garaldi E., Kannan R., Smith A., Springel V., Pakmor R., Vogelsberger M., Hernquist L., 2022, MNRAS, 512, 4909\n- Garcia F. A. B., Ricotti M., Sugimura K., Park J., 2023, MNRAS, 522, 2495 Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MNRAS, 448, 3248\n- Gennaro M., et al., 2018, ApJ, 863, 38\n- Gessey-Jones T., et al., 2022, MNRAS, 516, 841\n- Gnedin N. Y., Abel T., 2001, New A, 6, 437\n- Goldsmith P. F., Langer W. D., 1978, ApJ, 222, 881\n- Gong M., Ostriker E. C., Wolfire M. G., 2017, ApJ, 843, 38\n- Götberg Y., de Mink S. E., Groh J. H., 2017, A&A, 608, A11\n- Grudi'c M. Y., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Murray N., Kereš D., 2018, MNRAS, 475, 3511\n- Grudi'c M. Y., Guszejnov D., Hopkins P. F., O ff ner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 506, 2199\n- Gutcke T. A., Springel V., 2019, MNRAS, 482, 118\n- Gutcke T. A., Pakmor R., Naab T., Springel V., 2021, MNRAS, 501, 5597\n- Gutcke T. A., Pfrommer C., Bryan G. L., Pakmor R., Springel V., Naab T., 2022,\n- ApJ, 941, 120\n- Haardt F., Madau P., 2012, ApJ, 746, 125\n- Habing H. J., 1968, Bull. Astron. Inst. Netherlands, 19, 421\n- Haid S., Walch S., Seifried D., Wünsch R., Dinnbier F., Naab T., 2018, MNRAS,", '478, 4799': "- Harris C. R., et al., 2020, Nature, 585, 357\n- Hartmann L., Ballesteros-Paredes J., Heitsch F., 2012, MNRAS, 420, 1457\n- Hayashi K., Chiba M., Ishiyama T., 2020, ApJ, 904, 45\n- Hayward C. C., Hopkins P. F., 2017, MNRAS, 465, 1682\n- Heays A. N., Bosman A. D., van Dishoeck E. F., 2017, A&A, 602, A105\n- Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ, 591, 288\n- Hindmarsh A. C., Brown P. N., Grant K. E., Lee S. L., Serban R., Shumaker D. E., Woodward C. S., 2005, ACM Transactions on Mathematical Software (TOMS), 31, 363  \nHirai Y., Fujii M. S., Saitoh T. R., 2021, PASJ, 73, 1036  \n- Hislop J. M., Naab T., Steinwandel U. P., Lahén N., Irodotou D., Johansson P. H., Walch S., 2022, MNRAS, 509, 5938\n- Hollenbach D., McKee C. F., 1979, ApJS, 41, 555\n- Hollyhead K., Bastian N., Adamo A., Silva-Villa E., Dale J., Ryon J. E., Gazak Z., 2015, MNRAS, 449, 1106\n- Hopkins P. F., Quataert E., Murray N., 2011, MNRAS, 417, 950\n- Hopkins P. F., Quataert E., Murray N., 2012a, MNRAS, 421, 3488\n- Hopkins P. F., Quataert E., Murray N., 2012b, MNRAS, 421, 3522\n- Hopkins P. F., Narayanan D., Murray N., 2013, MNRAS, 432, 2647\n- Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581\n- Hopkins P. F., et al., 2018, MNRAS, 480, 800\n- Hopkins P. F., Grudi'c M. Y., Wetzel A., Kereš D., Faucher-Giguère C.-A., Ma X., Murray N., Butcher N., 2020, MNRAS, 491, 3702\n- Hopkins P. F., et al., 2023, MNRAS, 519, 3154\n- Hu C.-Y., 2019, MNRAS, 483, 3363\n- Hu C.-Y., Naab T., Walch S., Glover S. C. O., Clark P. C., 2016, MNRAS, 458, 3528\n- Hu C.-Y., Naab T., Glover S. C. O., Walch S., Clark P. C., 2017, MNRAS, 471, 2151\n- Hu C.-Y., Sternberg A., van Dishoeck E. F., 2021, ApJ, 920, 44\n- Hunter J. D., 2007, Computing in Science & Engineering, 9, 90\n- Indriolo N., McCall B. J., 2012, ApJ, 745, 91\n- Jaura O., Glover S. C. O., Klessen R. S., Paardekooper J. P., 2018, MNRAS, 475, 2822\n- Je ff reson S. M. R., Krumholz M. R., Fujimoto Y., Armillotta L., Keller B. W., Chevance M., Kruijssen J. M. D., 2021, MNRAS, 505, 3470\n- Jeon M., Besla G., Bromm V., 2017, ApJ, 848, 85\n- Jeon M., Bromm V., Besla G., Yoon J., Choi Y., 2021, MNRAS, 502, 1\n- Kannan R., Stinson G. S., Macciò A. V., Brook C., Weinmann S. M., Wadsley J., Couchman H. M. P., 2014, MNRAS, 437, 3529\n- Kannan R., Vogelsberger M., Marinacci F., McKinnon R., Pakmor R., Springel V., 2019, MNRAS, 485, 117\n- Kannan R., Marinacci F., Simpson C. M., Glover S. C. O., Hernquist L., 2020a, MNRAS, 491, 2088\n- Kannan R., Marinacci F., Vogelsberger M., Sales L. V., Torrey P., Springel V., Hernquist L., 2020b, MNRAS, 499, 5732\n- Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M., Hernquist L., 2022, MNRAS, 511, 4005\n- Kannan R., et al., 2023, MNRAS, 524, 2594\n- Karakas A. I., 2010, MNRAS, 403, 1413\n- Katz N., 1992, ApJ, 391, 502\n- Katz H., et al., 2022, MNRAS submitted, p. arXiv:2211.04626\n- Kennicutt Robert C. J., 1998, ApJ, 498, 541\n- Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99\n- Kim J., et al., 2021a, MNRAS, 504, 487\n- Kim J.-G., Ostriker E. C., Filippova N., 2021b, ApJ, 911, 128\n- Kim J.-G., Gong M., Kim C.-G., Ostriker E. C., 2023a, ApJS, 264, 10\n- Kim J., Jeon M., Choi Y., Richstein H., Sacchi E., Kallivayalil N., 2023b, ApJ, 959, 31\n- Kim C.-G., et al., 2024, ApJ in press, p. arXiv:2405.19227\n- Kroupa P., Weidner C., Pflamm-Altenburg J., Thies I., Dabringhausen J., Marks M., Maschberger T., 2013, in Oswalt T. D., Gilmore G., eds, , Vol. 5, Planets, Stars and Stellar Systems. Volume 5: Galactic Structure and Stellar Populations. p. 115, doi:10.1007 / 978-94-007-5612-0\\_4  \nKruijssen J. M. D., et al., 2019, Nature, 569, 519  \n- Krumholz M. R., McKee C. F., Bland-Hawthorn J., 2019, ARA&A, 57, 227\n- Lada C. J., Lada E. A., 2003, ARA&A, 41, 57\n- Lada C. J., Forbrich J., Lombardi M., Alves J. F., 2012, ApJ, 745, 190\n- Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., 2019, ApJ, 879, L18\n- Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., Steinwandel U. P., Moster B. P., 2020, ApJ, 891, 2  \nLahén N., et al., 2023, MNRAS, 522, 3092  \n- Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021a, ApJ, 914, 89\n- Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021b, ApJ, 914, 90\n- Lanz T., Hubeny I., 2003, ApJS, 146, 417\n- Larkin M. M., Gerasimov R., Burgasser A. J., 2023, AJ, 165, 2\n- Larson R. B., 2003, in De Buizer J. M., van der Bliek N. S., eds, Astronomical Society of the Pacific Conference Series Vol. 287, Galactic Star Formation Across the Stellar Mass Spectrum. pp 65-80 ( arXiv:astro-ph/0205466 ), doi:10.48550 / arXiv.astro-ph / 0205466\n- Lee T., Jeon M., Bromm V., 2024, MNRAS, 527, 1257\n- Leroy A. K., Walter F., Brinks E., Bigiel F., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2782\n- Levermore C. D., 1984, J. Quant. Spectr. Rad. Transf., 31, 149\n- Li H., Vogelsberger M., Marinacci F., Gnedin O. Y., 2019, MNRAS, 487, 364 Limongi M., Chie ffi A., 2018, ApJS, 237, 13\n- Lodders K., 2019, arXiv e-prints, p. arXiv:1912.00844\n- Maoz D., Mannucci F., Brandt T. D., 2012, MNRAS, 426, 3282\n- Marinacci F., Sales L. V., Vogelsberger M., Torrey P., Springel V., 2019, MNRAS, 489, 4233\n- Marks M., Kroupa P., Dabringhausen J., Pawlowski M. S., 2012, MNRAS, 422, 2246\n- Martin-Alvarez S., Sijacki D., Haehnelt M. G., Farcy M., Dubois Y., Belokurov\n- V., Rosdahl J., Lopez-Rodriguez E., 2023, MNRAS, 525, 3806  \nMateo M. L., 1998, ARA&A, 36, 435  \n- McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148\n- McKinnon R., Vogelsberger M., Torrey P., Marinacci F., Kannan R., 2018, MNRAS, 478, 2851\n- McKinnon R., Kannan R., Vogelsberger M., O'Neil S., Torrey P., Li H., 2021, MNRAS, 502, 1344\n- Meynet G., Maeder A., 2005, A&A, 429, 581\n- Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121\n- Nomoto K., Iwamoto K., Nakasato N., Thielemann F. K., Brachwitz F., Tsujimoto T., Kubo Y., Kishimoto N., 1997, Nucl. Phys. A, 621, 467\n- Ostriker E. C., Kim C.-G., 2022, ApJ, 936, 137\n- Ostriker E. C., Shetty R., 2011, ApJ, 731, 41\n- Ostriker E. C., McKee C. F., Leroy A. K., 2010, ApJ, 721, 975\n- Pakmor R., Bauer A., Springel V., 2011, MNRAS, 418, 1392\n- Pakmor R., Springel V., Bauer A., Mocz P., Munoz D. J., Ohlmann S. T., Schaal\n- K., Zhu C., 2016, MNRAS, 455, 1134\n- Parmentier G., Kau ff mann J., Pillai T., Menten K. M., 2011, MNRAS, 416, 783\n- Peter T., Klessen R. S., Kanschat G., Glover S. C. O., Bastian P., 2023, MNRAS, 519, 4263\n- Petkova M., Springel V., 2009, MNRAS, 396, 1383\n- Pillepich A., et al., 2018, MNRAS, 473, 4077\n- Ploeckinger S., Schaye J., 2020, MNRAS, 497, 4857\n- Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431 Portinari L., Chiosi C., Bressan A., 1998, A&A, 334, 505\n- Prgomet M., Rey M. P., Andersson E. P., Segovia Otero A., Agertz O., Renaud F., Pontzen A., Read J. I., 2022, MNRAS, 513, 2326\n- Rahmati A., Pawlik A. H., Raiˇcevi'c M., Schaye J., 2013, MNRAS, 430, 2427 Rathjen T.-E., et al., 2021, MNRAS, 504, 1039\n- Rémy-Ruyer A., et al., 2014, A&A, 563, A31\n- Renaud F., Romeo A. B., Agertz O., 2021, MNRAS, 508, 352\n- Revaz Y., Arnaudon A., Nichols M., Bonvin V., Jablonka P., 2016, A&A, 588, A21\n- Richings A. J., Schaye J., 2016, MNRAS, 458, 270\n- Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188\n- Rosdahl J., Schaye J., Teyssier R., Agertz O., 2015, MNRAS, 451, 34  \nRosen A., Bregman J. N., 1995, ApJ, 440, 634  \n- Roychowdhury S., Huang M.-L., Kau ff mann G., Wang J., Chengalur J. N., 2015, MNRAS, 449, 3700\n- Schaerer D., 2002, A&A, 382, 28\n- Schaye J., 2004, ApJ, 609, 667\n- Schaye J., et al., 2015, MNRAS, 446, 521\n- Schinnerer E., Leroy A. K., 2024, ARA&A in press, p. arXiv:2403.19843\n- Schmidt M., 1959, ApJ, 129, 243\n- Secunda A., Cen R., Kimm T., Götberg Y., de Mink S. E., 2020, ApJ, 901, 72\n- Semenov V. A., Kravtsov A. V., Gnedin N. Y., 2017, ApJ, 845, 133\n- Shi Y., Helou G., Yan L., Armus L., Wu Y., Papovich C., Stierwalt S., 2011, ApJ, 733, 87\n- Shi Y., et al., 2018, ApJ, 853, 149\n- Skinner M. A., Ostriker E. C., 2013, ApJS, 206, 21\n- Smith N., 2014, ARA&A, 52, 487\n- Smith M. C., 2021, MNRAS, 502, 5417\n- Smith B. D., et al., 2017, MNRAS, 466, 2217\n- Smith A., Kannan R., Tsang B. T. H., Vogelsberger M., Pakmor R., 2020, ApJ, 905, 27\n- Smith M. C., Bryan G. L., Somerville R. S., Hu C.-Y., Teyssier R., Burkhart B., Hernquist L., 2021, MNRAS, 506, 3882\n- Smith A., Kannan R., Garaldi E., Vogelsberger M., Pakmor R., Springel V., Hernquist L., 2022, MNRAS, 512, 3243\n- Springel V., 2010, MNRAS, 401, 791\n- Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776\n- Steinwandel U. P., Goldberg J. A., 2023, ApJ submitted, p. arXiv:2310.11495\n- Steinwandel U. P., Moster B. P., Naab T., Hu C.-Y., Walch S., 2020, MNRAS, 495, 1035\n- Steinwandel U. P., Bryan G. L., Somerville R. S., Hayward C. C., Burkhart B., 2023, MNRAS, 526, 1408\n- Steinwandel U. P., Rennehan D., Orr M. E., Fielding D. B., Kim C.-G., 2024, ApJ submitted, p. arXiv:2407.14599  \nSu K.-Y., et al., 2018, MNRAS, 480, 1666  \n- Sugimura K., Ricotti M., Park J., Garcia F. A. B., Yajima H., 2024, ApJ, 970, 14 Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H. T., 2016, ApJ, 821, 38\n- Sun J., et al., 2023, ApJ, 945, L19", 'Appendix A: Treatments of cooling channels': 'In this appendix, we provide a detailed description of the various cooling channels in our model.', 'Appendix A.1: Primordial cooling': "The term 'primordial cooling' incorporates photoionization, recombination, free-free (bremsstrahlung), bound-bound and bound-free collisions of the H and He species tracked in the nonequilibrium network. In this work, we use the primordial thermochemistry and cooling equations described by equations (49)-(51) of Kannan et al. (2019) and equations (2)-(5) of Kannan et al. (2020b). We adopt the rate coe ffi cients accounting for cooling / heating of atomic hydrogen and helium from Cen (1992). We model the cooling due to molecular hydrogen by  \nn 2 Λ H2 = Λ H i H2 , n → 0( T ) n H i n H2 + Λ H2H2 , n → 0( T ) n H2 n H2 , (A.1)  \nwhere Λ H i H2 , n → 0 and Λ H2H2 , n → 0 are the H i -H2 and H2-H2 collisional cooling coe ffi cients in the low-density limit ( n → 0) from Hollenbach & McKee (1979). We temporarily ignore UV pumping heating and H2 formation heating as they are subdominant in our current framework.", 'Appendix A.2: Metal collisional excitation lines in hot gas': 'The metal-line cooling in ≳ 10 5 K hot gas is modelled assuming ionization equilibrium in a UVB radiation field given by Faucher-Giguère et al. (2009) with a density-dependent selfshielding factor from Rahmati et al. (2013). Given the density and temperature of a gas cell, we look up this cooling rate from a pre-calculated cloudy table (see Vogelsberger et al. 2013 for more details).', 'Appendix A.3: Nebular lines in photoionized gas': "In photoionized gas with sub-solar abundance ( Z > 0 . 1Z ⊙ ), the collisionally excited forbidden lines of metal ions become important coolants. We adopt the fitting formula given by equation (47) of Kim et al. (2023a) for the nebular lines except those from C + assuming a fixed ionization state such that 80% (20%) of O, N, and Ne are singly (doubly) ionized and 50% (50%) of S is singly (doubly) ionized:  \nn 2 Λ neb = 3 . 86 × 10 -23 erg cm -3 s -1 × Z ' g nen H ii e -3 . 86 / T 4 f neb( T ) T 1 / 2 4 (1 + 0 . 12 n 0 . 38 -0 . 12 ln T 4 e , 2 ) , (A.2)  \nwhere Z ' g is the gas metallicity normalized to the solar value, T 4 = T / 10 4 K, ne , 2 = ne / 10 2 cm -3 , and log 10 f neb( T ) = P 6 i = 0 ai (ln T 4) i with { a 0 , a 1 , a 2 , a 3 , a 4 , a 5 } = { 0 . 692 , -0 . 586 , 0 . 816 , -0 . 505 , 0 . 118 , 7 . 66 × 10 -3 , -5 . 08 × 10 -3 } . This function captures the temperature dependence of the collisional excitation rate of the dominant nebular coolant [O ii ]3729 Å line and the critical densities of multiple lines for their collisional deexcitation in dense gas (for more details, see Kim et al. 2023a).", 'Appendix A.4: Metal fine structure lines cooling in warm gas': 'The main metal coolants of warm atomic gas are fine-structure transitions of C + and O with minor contributions from C and  \nFig. A.1. Cooling functions in predominantly neutral gas with n H = 10 cm -3 , xe = 2 . 7 × 10 -3 , G 0 = 1 for di ff erent metallicities of Z = 1 Z ⊙ (top panel), 0 . 1 Z ⊙ (middle panel), and 0 . 01 Z ⊙ (bottom panel) obtained with the FIRE-2 (green), FIRE-3 (blue), and RIGEL (black) model.  \n<!-- image -->  \nT  \n[K]  \nother species (e.g. Wolfire et al. 1995, 2003). As noticed by Kim et al. (2023a), the low-temperature ( < 10 4 K) metal cooling models widely used by the galaxy formation community have significant discrepancies with those developed with detailed atomic / molecular physics by the ISM community. To accurately model such low-temperature metal cooling, we calculate the cooling rates of the predominant C / O fine structure lines including C ii 158 µ m, C i ∗ 610 µ m, C i ∗ 370 µ m, O i ∗ 63 µ m, and O i ∗ 146 µ m lines. To obtain the cooling rate of a given emission line, we calculate the equilibrium abundances of C, C + , CO, O, and O + under the radiation field and nonequilibrium H and He abundances tracked by arepo -rt (see Section 2.2), then we solve the two-/ three-level population of C / O atoms / ions with the rate coe ffi cients given by Gong et al. (2017). For example, the cooling rate for the C ii 158 µ mtwo-level system is  \nn 2 Λ C ii = n C ii A ul E ul njk lu njk lu + njk ul + A ul , (A.3)  \nwhere A ul is the Einstein A coe ffi cient, E ul is the energy di ff erence, k ul is the collisional rate coe ffi cien and k lu = g u / g l k ul exp ( -E ul / k B T ), and nj is the number density of collisional partners including H2, H i , and free electron. The cooling rates for the three-level systems are calculated in the same way and the solution for the three-level population can be found in Section 17.5 of Draine (2011).  \nIn Fig. A.1, we compare the cooling functions in predominantly neutral gas with n H = 10 cm -3 , xe = 2 . 7 × 10 -3 , G 0 = 1  \nobtained from our RIGEL model with those from the FIRE-2 (Hopkins et al. 2018) and FIRE-3 (Hopkins et al. 2023) models as an example to show the di ff erence. The cooling curves in different metallicities of 1 Z ⊙ , 0 . 1 Z ⊙ , and 0 . 01 Z ⊙ are shown in the panels from top to bottom. Both RIGEL and FIRE-3 show a strong dependence on metallicity in their cooling rates, while this dependence in the FIRE-2 model is notably weak. Consequently, the FIRE-2 model can significantly overestimate the cooling rate in metal-poor environments. This weak metallicity dependency in FIRE-2 can be attributed to the inclusion of terms such as [1 + ( Z / Z ⊙ )] / (1 + 0 . 001 43 n H) in their fitting function.  \nCompared to FIRE-3, RIGEL demonstrates overall higher cooling rates. This variation occurs due to two main factors. First, RIGEL incorporates the cooling e ff ects from the O i metal fine structure lines alongside the C i and C ii lines present in the FIRE-3 model. Second, the species abundance and collisional rate coe ffi cients from all these lines display intricate dependencies on temperature, the density of the collision partners, and the radiation field, while FIRE-3 simplified these dependencies to obtain a simple fitting function.', 'Appendix A.5: CO cooling in cold gas': 'In dense ( n ≳ 10 3 cm -3 ) molecular gas, the rotational lines of CO and other molecules become the most important coolants. We model the cooling by CO rotational lines following the approach outlined by Whitworth & Ja ff a (2018),  \nn 2 Λ CO = n CO x H2 GLYPH<16> Λ -1 /β CO,LO + Λ -1 /β CO,HI GLYPH<17> -β , (A.4)  \nwhere Λ CO,LO, Λ CO,HI, and β are given by  \n3 2  \nΛ CO,LO = 2 . 16 × 10 -27 erg s -1 GLYPH<18> n H2 cm -3 GLYPH<19> GLYPH<18> T K GLYPH<19> / , Λ CO,HI = 2 . 21 × 10 -28 erg s -1 × " ( x CO / x H2) km s -1 pc -1 |∇ · v | # -1 GLYPH<18> n H2 cm -3 GLYPH<19> -1 GLYPH<18> T K GLYPH<19> 4 , β = 1 . 23 GLYPH<18> n H2 cm -3 GLYPH<19> 0 . 0533 GLYPH<18> T K GLYPH<19> 0 . 164 . (A.5)  \nwhich reproduces the detailed numerical results of Goldsmith & Langer (1978) for T ∈ [10 , 100] K, n H2 ∈ h 3 × 10 2 , 3 × 10 5 i cm -3 , and |∇ · v | [ x CO / x H2 / (3 × 10 -4 )] -1 ∈ h 10 , 10 6 i kms -1 pc -1 .', 'Appendix A.6: Dust cooling in cold gas': "In even denser molecular gas ( n ≳ 10 6 cm -3 ) where molecules reach their critical densities, the dust-gas energy exchange can be an important coolant. The dust cooling function is given by Burke & Hollenbach (1983) and rearranged by Kannan et al. (2020b):  \nn 2 Λ gd = 1 . 356 × 10 -33 erg cm -3 s -1 × GLYPH<18> Z d 0 . 01 GLYPH<19> 0 . 1 µ m a ! r T K GLYPH<18> T -T d K GLYPH<19> GLYPH<18> n H cm -3 GLYPH<19> 2 . (A.6)  \nwhere T d is the dust temperature which is self-consistently calculated considering the dust-radiation coupling. The size of the dust grain a is fixed at 0 . 1 µ m. We calculate T d by solving the  \nArticle number, page 28 of 30  \ninstantaneous equilibrium equation of Λ gd +Γ dr = 0 using Newton's root-finding method. The dust heating rate due to the coupling with IR radiation field is given by  \nn 2 Γ dr = κ P ρ (4 σ T 4 d -cE IR) , (A.7)  \nwhere κ P is the Planck mean opacity, σ is the Stefan-Boltzmann constant, and E IR is the energy density of IR photons.", 'Appendix A.7: Photoelectric heating': "FUV photons absorbed by small dust grains and polycyclic aromatic hydrocarbons (PAHs) can excite electrons to escape and then heat the ISM, which is referred to as photoelectric heating. As in Kannan et al. (2020b), we adopt the photoelectric heating rate given by Wolfire et al. (2003)  \nn 2 Γ PE = 1 . 3 × 10 -24 erg s -1 cm -3 GLYPH<18> Z d 0 . 01 GLYPH<19> ϵ ff G 0 n H , (A.8)  \nwhere G 0 is the FUV radiation energy density relative to Habing (1968)'s estimate ( u 5.8-13.6 eV = 5 . 29 × 10 -14 erg cm -3 ) and  \nϵ ff = 0 . 0487 1 + 0 . 004[ G 0 √ T / (0 . 5 ne )] 0 . 73 + 0 . 0365( T / 10 4 ) 0 . 7 1 + 2 × 10 -4 × G 0 √ T / (0 . 5 ne ) . (A.9)  \nHere, the unit of temperature T is [K], and the unit of electron number density ne is [cm -3 ].", 'Appendix A.8: Cosmic ray heating': 'We assume constant ionization and heating rates with 0.01 times of the MW value, where the MW values (Indriolo & McCall 2012) are  \nζ cr,H2 = 3 . 5 × 10 -16 s -1 ζ cr,H i = 1 . 78 × 10 -16 s -1 ζ cr,He i = 1 . 1 ζ cr,H i ζ cr,C i = 3 . 85 ζ cr,H i n 2 Γ cr , j = 1 . 6022 × 10 -11 erg · ζ cr , jnj , (A.10)  \nwhere j ∈ H2 , H i , He i .', 'Appendix B: Tables of fitting parameters': 'In Table B.1, Table B.2, and Table B.3, we listed the fitting parameters x 0 and bi ( i = 0 , 1 , 2 , 3 , 4) the for main-sequence lifetime (in [yr]), UV photon production rates in di ff erent bands (in [s -1 ]), and wind mass-loss rate (in [M ⊙ yr -1 ]) and velocity (in [kms -1 ]), defined in Eq. (8), respectively. The machine-readable tables with full precision are downloadable online 5 .  \nTable B.1. Fit parameters for main-sequence lifetime of stars.  \nNotes. We adopt the Portinari et al. (1998) model for Pop. II stars and the S02 model for Pop. III ( Z ≤ 10 -8 Z ⊙ ) stars. We use the fit formula of Eq. 8 fixing b 4 = 0. The unit of main-sequence lifetime is [yr].  \nTable B.2. Fit parameters for photon production rates in di ff erent bands (see Table 1).  \nNotes. Opt. (1 -5 . 8 eV), FUV (5 . 8 -11 . 2 eV), LW (11 . 2 -13 . 6 eV), EUV1 (13 . 6 -24 . 6 eV), EUV2 (24 . 6 -54 . 4 eV), EUV3 ( > 54 . 4 eV). The unit of photon production rate is [s -1 ].  \nTable B.3. Fit parameters for wind mass-loss rate and velocity.  \nNotes. The wind properties is calculated by Vink et al. (2001); Vink & Sander (2021) with the ZAMS stellar properties from T20 including the extrapolated model for Z = Z ⊙ . We use the fit formula of Eq. (8). The units of mass-loss rate and wind velocity are [M ⊙ yr -1 ] and [km s -1 ], respectively.', 'Appendix C: Resolution of SN explosions': 'To quantitatively assess the ability of our simulations to resolve the ST phase of SN explosion, we record the mass and gas density of the cell where SN energy is ejected. In Fig. C.1, we present a histogram of the ratio between the mass of the cell where the SN has injected explosion energy and the shell formation mass described by equation (12).  \nAs can be seen, the fiducial RIGEL model and the noRT model at 1 M ⊙ resolution show a good ability to resolve the SNe. For the 10 M ⊙ low-resolution simulations, the presence of radiative feedback ensures that all SN explosions occur in lowdensity gas (Section 3.4.4), where the 10 M ⊙ resolution is su ffi -cient to resolve all SNe. In the 0 . 02Z ⊙ / noRT / low simulation, a small fraction of SNe is marginally resolved, while they should have minor e ff ects on the results.  \nFig. C.1. Cumulative distribution function of the ratio between the mass of the cell where the SNe injected explosion energy and the total sweptup mass at end of the energy conserving Sedov-Taylor phase (shell formation) given by equation (12). The red and green curves are the 1 M ⊙ resolution results obtained with the fiducial RIGEL model, while the blue curve is obtained without radiative feedback (see Section 3 for the naming conventions). The dotted curves are the results of the 10 M ⊙ resolution simulations. The vertical dashed line denotes the STresolving criterion of m gas,SN / M shell < 1 / 27 by Kim & Ostriker (2015).  \n<!-- image -->'}