© The Authors 2024

1. Introduction

Star formation on galactic scales is remarkably inefficient (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 inefficiency is particularly pronounced in dwarf galaxies (M_⋆ < 10⁹ 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 efficiency (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 deposit 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, 2012a).

Stellar 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ć et al. 2018; Li et al. 2019; Kim et al. 2021a). 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 that of SN explosions (Agertz et al. 2013; Kannan et al. 2014; Geen et al. 2015; Hopkins et al. 2020; Jeffreson 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. 2021b). 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 effects 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).

Over the past few decades, cosmological simulations have been able to reproduce the inefficiency 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).

Recent advances in computational capabilities have enabled simulations of the cold (≲10⁴ K) dense (≳1 cm⁻³) 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. 2020a). 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³ 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 effective photoionization feedback, and an incorrect baryonic cycle within these smaller galaxies (Su et al. 2018; Applebaum et al. 2020; Smith 2021).

For 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; Applebaum 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⁷ 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. 2012b, 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 effects on the formation and population of star clusters.

Although 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 effective 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 sufficiently accurate in the wide variety of realistic environments encountered in simulations (Rosdahl et al. 2015; Fujimoto et al. 2019; Kannan et al. 2020b). The lack of long-range radiative effects in these effective 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 affordably 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 efficiency and accuracy of AREPO-RT (Kannan et al. 2019; Deng et al. 2024), now on-the-fly RT simulations are already obtainable for intermediate-sized dwarf galaxies (e.g., M_gas = 4 × 10⁷ M_⊙) with a high resolution of ∼1 M_⊙.

In 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 of boosted radiation pressure due to increased opacity in their atmospheres. These metallicity-dependent factors can introduce orders-of-magnitude differences in the feedback budget of stars formed in either enriched or pristine environments, which must be reflected in cosmological simulations.

Finally, 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⁴ 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 efficient 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.

In this paper, we introduce the Realistic ISM modeling in Galaxy Evolution and Lifecycles (RIGEL) model, a self-consistent model for simulating stellar feedback and the multiphase ISM in dwarf galaxies. Our model brings together on-the-fly 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 effects of turning RT on or off. We discuss the implications and caveats of our model in Section 4. Lastly, we give a summary in Section 5.

2. The RIGEL model

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¹). 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 octree 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.

The 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, $\stackrel{\sim }{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 second-order 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. (2020a), with an update to the low-temperature (≲10⁴ 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⁺, H₂, He, He⁺, and He²⁺ – through AREPO-RT. We refer the readers to Section 2.1 in Kannan et al. (2020a) 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).

We also followed the enrichment and evolution of seven individual metal species including C, N, O, Ne, Mg, Si, and Fe². The metals released by the stars diffused and mixed numerically 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. (2020a), which treats metal cooling as tabulated cooling rates from CLOUDY, we modeled high-temperature (≳10⁵ K) and low-temperature (≲10⁴ 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.5Z, 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 ${\zeta }_{\mathrm{cr},\mathrm{H}\mathrm{I}}=1.78\times {10}^{-18}{\mathrm{s}}^{-1}$, and ζ_cr is used for ${\zeta }_{\mathrm{cr},\mathrm{H}\mathrm{I}}$ if not specified.

To 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

$$\begin{array}{c}\hfill \frac{x_{{\mathrm{C}}^{+}}}{x_{\mathrm{C},\mathrm{tot}}}=\frac{{\zeta }_{\mathrm{pi},\mathrm{C}}+{\zeta }_{\mathrm{cr},\mathrm{C}}}{{\zeta }_{\mathrm{pi},\mathrm{C}}+{\zeta }_{\mathrm{cr},\mathrm{C}}+{\alpha }_{\mathrm{gr},{\mathrm{C}}^{+}}n+{\alpha }_{\mathrm{rr}+\mathrm{dr},{\mathrm{C}}^{+}}{n}_e+k_{{\mathrm{C}}^{+}-{\mathrm{H}}_2}{n}_{\mathrm{H}2}},\end{array}$$(1)

where we have assumed that $x_{\mathrm{C}}=x_{\mathrm{C},\mathrm{tot}}-x_{{\mathrm{C}}^{+}}$, since x_CO is negligible in C/C⁺ transition regions.

In our model, C is ionized by photoionization and CRs. The photoionization rate, ${\zeta }_{\mathrm{pi},\mathrm{C}}={\zeta }_{\mathrm{pi},\mathrm{C}}^{\mathrm{Draine}}{\chi }_{\mathrm{LW}}{f}_{\mathrm{s},\mathrm{C}}+520\times 2x_{{\mathrm{H}}_2{\zeta }_{\mathrm{cr}}}$, where χ_LW is the LW band radiation intensity (tracked by AREPO-RT) normalized by the Draine (1978) value of 1.3 × 10⁻¹⁴ erg cm⁻³, ${\zeta }_{\mathrm{pi},\mathrm{C}}^{\mathrm{Draine}}=3.43\times {10}^{-10}$ s⁻¹ 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 coefficients for the formation channels of C, including gas phase and grain surface recombination (${\alpha }_{\mathrm{rr}+\mathrm{dr},{\mathrm{C}}^{+}}$ and ${\alpha }_{{\mathrm{gr},\mathrm{C}}^{+}}$) and the reaction of ${\mathrm{C}}^{+}+{\mathrm{H}}_2\to {\mathrm{CH}}_2^{+}$ ($k_{{\mathrm{C}}^{+}-{\mathrm{H}}_2}$), were taken from Gong et al. (2017).

The O⁺ abundance was determined by assuming $x_{{\mathrm{O}}^{+}}/x_{\mathrm{O},\mathrm{tot}}=x_{\mathrm{H}\mathrm{I}\mathrm{I}}$ 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_{{\mathrm{C}}^{+}}$ and $x_{\mathrm{O}}^{+}$ were determined, we determined the CO fraction x_CO by

$$\begin{array}{c}\hfill \frac{x_{\mathrm{CO}}}{x_{\mathrm{C},\mathrm{tot}}}=\frac{2x_{{\mathrm{H}}_2}[1-max(x_{{\mathrm{C}}^{+}}/x_{\mathrm{C},\mathrm{tot}},x_{{\mathrm{O}}^{+}}/x_{\mathrm{O},\mathrm{tot}})]}{1+{(n_{\mathrm{CO},\mathrm{crit}}/n)}^2},\end{array}$$(2)

where 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),

$$\begin{array}{c}\hfill n_{\mathrm{CO},\mathrm{crit}}={(4\times {10}^{3}Z{\prime }_{\mathrm{d}}{\zeta }_{\mathrm{cr},16}^{-2})}^{{\chi }_{\mathrm{LW}}^{1/3}}\left(\frac{50{\zeta }_{\mathrm{cr},16}}{{Z\prime }_{\mathrm{d}}^{1.4}}\right){\mathrm{cm}}^{-3},\end{array}$$(3)

where ζ_cr, 16 = ζ_cr/(10⁻¹⁶ s⁻¹) and Z′_d is the dust metallicity normalized by the solar value.

Based on the chemical abundance of the main ISM coolants, we can calculate the cooling rate of different 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⁵ 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

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_{\mathrm{net}}=& {\mathrm{\Lambda }}_{\mathrm{pri}}(n_j,N_{\gamma }^i,T)+\frac{Z}{{\mathrm{Z}}_{\odot }}{\mathrm{\Lambda }}_{\mathrm{Z},\mathrm{CIE}}(T,\rho )\hfill \\ \hfill & +{\mathrm{\Lambda }}_{\mathrm{Z},\mathrm{C}/\mathrm{O}}(n_j,n_{\mathrm{C}/\mathrm{O}},N_{\gamma }^{\mathrm{LW}},T)\hfill \\ \hfill & +{\mathrm{\Lambda }}_{\mathrm{Z},\mathrm{neb}}(n_{\text{H}\text{II}},n_e,T,Z_{\mathrm{g}})+{\mathrm{\Lambda }}_{\mathrm{dust}}(\rho ,T,D,N_{\gamma }^{\mathrm{IR}},Z_{\mathrm{d}})\hfill \\ \hfill & -{\mathrm{\Gamma }}_{\mathrm{PE}}(D,T,N_{\gamma }^{\mathrm{FUV}},Z_{\mathrm{g}})-{\mathrm{\Gamma }}_{\mathrm{cr}}(n_j),\hfill \end{array}$$(4)

where n_j is the number density of all the primordial species ($j\in \{{\mathrm{H}}_2,\mathrm{H}\mathrm{I},\mathrm{H}\mathrm{I}\mathrm{I},\mathrm{He}\mathrm{I},\mathrm{He}\mathrm{I}\mathrm{I},\mathrm{He}\mathrm{I}\mathrm{I}\mathrm{I}\}$), N_γⁱ 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.

Closely 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|² + |∇×v|² < 2Gρ, Hopkins et al. 2013), and marginally Jeans-resolved ($m_{\mathrm{gas}}<f_{\mathrm{J},\mathrm{n}}^3{M}_{\mathrm{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 𝒫_SF = Δt/τ_ff, where τ_ff = (3π/32Gρ)^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_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 different definitions of M_J and selections of f_J, s in the literature. We followed the numerical experiments by Grudić et al. (2021), which suggest that f_J, s = 0.5 is acceptable with the Jeans mass defined as $M_{\mathrm{J}}={\pi }^{1.5}{c}_{\mathrm{s}}^3{G}^{-1.5}{\rho }^{-0.5}$. The parameters used for the simulations in this work are presented in Section 3.1.

Newly 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 differs from the dynamical mass of the star particles, m_⋆, and we determined the value of M_fb using the method outlined below.

Given the IMF, Φ(M, Z)≡dN/dM, 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, 𝒫_host, of hosting a massive star,

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{host}}=m_{\star }{\displaystyle {\int }_{M_{\min ,\mathrm{SN}}}^{M_{\max ,\mathrm{IMF}}}}\mathrm{\Phi }(M,Z)\mathrm{d}M/{\displaystyle {\int }_{M_{\min ,\mathrm{IMF}}}^{M_{\max ,\mathrm{IMF}}}}M\mathrm{\Phi }(M,Z)\mathrm{d}M,\end{array}$$(5)

where 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 high-resolution simulations with m_⋆ ≲ 10 M_⊙, 𝒫_host < 1 always holds for the considered range of IMF variations.

If 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⁻¹(y), where y is a random variable with the value of [0, 1] and F⁻¹ is the inverse function of

$$\begin{array}{c}\hfill F(x)\equiv {\displaystyle {\int }_{M_{\min ,\mathrm{SN}}}^x}\mathrm{\Phi }(M,Z)\mathrm{d}M/\left[{\displaystyle {\int }_{M_{\min ,\mathrm{SN}}}^{M_{\max ,\mathrm{IMF}}}}\mathrm{\Phi }(M,Z)\mathrm{d}M\right].\end{array}$$(6)

If the IMF has a constant slope α (e.g., α = −2.3) at the M > M_min, SN range, F⁻¹(y) has a simple form of

$$\begin{array}{c}\hfill M_{\mathrm{fb}}={[(M_{\max }^{\alpha +1}-M_{\min }^{\alpha +1})y+M_{\min }^{\alpha +1}]}^{1/(\alpha +1)}.\end{array}$$(7)

For 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., metallicity-dependent IMF, see Kroupa et al. 2013 for a review) is straightforward under this framework and will be explored in future works.

We 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 of 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 different literature and convert them into functions of M and Z. We found that the metallicity dependence of the ZAMS stellar radius and effective temperature can be well approximated with a power law in the mass range of 8 − 120 M_⊙ for $Z\sim {10}^{-6}-0.01{\mathrm{Z}}_{\odot }$, 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).

For the convenience of implementation, we used a polynomial formula to fit the stellar properties, 𝒜(M, Z′), where 𝒜 can be the photon production rate in [s⁻¹], mass-loss rate in $[{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}]$, wind velocity in [km s⁻¹], 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₁₀(Z′) space. In practice, the interpolation over log₁₀(Z′) starts from a nonzero small metallicity Z₀, and the Z = Z₀ model is applied to all cases with Z < Z₀. We adopted Z₀ = 10⁻⁸ Z_⊙, below which stellar evolution is effectively metal-free³. The polynomial fitting formula has the form

$$\begin{array}{cc}& {log}_{10}\mathcal{A}=\{\begin{array}{cc}{b}_0+b_{1}x+b_2{x}^2+b_3{x}^3+b_4{x}^4,\hfill & x\le x_0,\hfill \\ {\mathcal{A}}_0+\frac{\mathrm{d}{log}_{10}\mathcal{A}}{\mathrm{d}x}{|}_{x=x_0}(x-x_0),\hfill & x>x_0,\hfill \end{array}\hfill \\ \hfill & {\mathcal{A}}_0=b_0+b_1{x}_0+b_2{x}_0^2+b_3{x}_0^3+b_4{x}_0^4,\hfill \\ \hfill & \frac{\mathrm{d}{log}_{10}\mathcal{A}}{\mathrm{d}x}{|}_{x=x_0}=(b_1+2b_2{x}_0+3b_3{x}_0^2+4b_4{x}_0^3)\hfill \end{array}$$(8)

for x ≡ log₁₀(M/M_⊙), where x₀ is a parameter chosen to optimize the fit quality at the high-mass end (which is better described by a linear relation between log₁₀𝒜 and log M). The relevant fit parameters, b_i (i = 0, 1, 2, 3, 4) and x₀, are shown in Appendix B (see Tables B.1, B.2, and B.3). In most cases, we fixed b₄ = 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 = 0.0004 − 0.05 combined with the function of metal-free Pop. III stars ($Z\le {10}^{-8}{\mathrm{Z}}_{\odot }$) 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.

Fig. 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}{\mathrm{Z}}_{\odot }$) 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}{\mathrm{Z}}_{\odot }$) to high (1 Z⊙). We shifted the curves of the Opt., FUV, LW, and EUV1 bands by multiplying constants of 104, 103, 102, and 10, respectively, to avoid overlapping different bins. Right: Similar to left but for the wind mass-loss rate (${\dot{m}}_{\mathrm{w}}$) and wind velocity (vw). The color of the curves from light to dark indicates the metal abundance from low (${10}^{-6}{\mathrm{Z}}_{\odot }$) to high (1 Z⊙). The main-sequence stellar winds from extremely metal-poor stars with $Z<{10}^{-6}{\mathrm{Z}}_{\odot }$ are negligible and not considered in our model.

During 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.

Once 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 & Chieffi (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⁻⁵, 3.11 × 10⁻⁴, 3.11 × 10⁻³, 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 intermediate-mass 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.

In 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 photodissociation, photoionization, and photoheating of H₂, 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.

The stellar spectra were discretized into seven bins (IR, Opt., FUV, LW, EUV1, EUV2, and EUV3) corresponding to our radiation bins. The photon production rate, Q_i, 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. 2020a for details). Similar to Kannan et al. (2020a), 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 Q_iΔt band i photons with mean energy e_i. The direction of the photon flux, F_γⁱ, was set to be radially outward from the star particle with a magnitude of $|{\mathit{F}}_{\gamma }^i|=\stackrel{\sim }{c}{e}_i{Q}_i\mathrm{\Delta }t/V_{\mathrm{cell}}$, where V_cell is the volume of the gas cell receiving the photons.

The photoionization and photon heating rates of species “j” are given by

$$\begin{array}{c}\hfill {\zeta }_{\mathrm{pi},j}=\stackrel{\sim }{c}{\displaystyle \sum _i}{\overline{\sigma }}_{\mathit{ij}}{N}_{\gamma }^i\end{array}$$(9)

and

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{pi},j}=\stackrel{\sim }{c}{\displaystyle \sum _i}{N}_{\gamma }^i{\overline{\sigma }}_{\mathit{ij}}{\mathfrak{h}}_{\mathit{ij}},\end{array}$$(10)

where ${\overline{\sigma }}_{\mathit{ij}}$ and 𝔥_ij are the mean ionization cross-section and photoelectron energy of photons in bin “i” interacting with species “j” (Kannan et al. 2019, 2020a).

Radiation 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

$$\begin{array}{c}\hfill {\frac{\partial \rho v}{\partial t}|}_{\mathrm{rad}.}=\frac{1}{c}{\displaystyle \sum _i}{\mathit{F}}_{\gamma }^i({\displaystyle \sum _j}{n}_j{\overline{\sigma }}_{\mathit{ij}}{p}_{\mathit{ij}}+{\kappa }_i\rho e_i),\end{array}$$(11)

where p_ij is the mean photon momentum and κ_i is the opacity due to dust (see Table 1 in Kannan et al. 2020a 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).

In general, the photon production rate, Q_i, 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 dust-gas coupling (Kannan et al. 2019, 2020a). For the six photon energy bands above 1 eV, the stellar evolution models⁴ used to construct the metallicity grids are described as follows.

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⁻⁸, 10⁻⁶, 10⁻⁵, 10⁻⁴, 0.01, and 1 Z_⊙ in T20 (EXT), assuming black-body spectra for simplicity.

2. For the LW band, the metallicity dependence is unclear in the literature. For simplicity, we considered three models at Z = 10⁻⁸, 10⁻⁴, 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}{\mathrm{Z}}_{\odot }$ ZAMS properties in 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⁻⁸ Z_⊙.

3. For the EUV1 and EUV2 bands that dominate the ionization feedback from massive stars, we considered four models at Z = 10⁻⁸, 0.001, 0.0282, and 1 Z_⊙. We used the Z = 0 model in S02 (see their Table 6) for Z ≤ 10⁻⁸ 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⁻⁴ 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.

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⁻⁸, 10⁻⁶, 10⁻⁵, 10⁻⁴, 0.01, and 1 Z_⊙.

We summarize the properties of the seven radiation bands in Table 1. In the middle panel of Fig. 1, we present the photon production rate, Q_i, 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

We modeled the main-sequence winds of massive stars by injecting the mass, ${\dot{m}}_{\mathrm{w}}\mathrm{\Delta }t$, and momentum, $\sqrt{2L_{\mathrm{w}}{\dot{m}}_{\mathrm{w}}}\mathrm{\Delta }t$, fluxes from a star to its nearest 32 neighboring gas cells weighted by their solid angle opening to the star (see Li et al. 2019 for details), where ${\dot{m}}_{\mathrm{w}}$ is the wind mass-loss rate and $L_{\mathrm{w}}=0.5{\dot{m}}_{\mathrm{w}}{v}^2$ is the wind luminosity or power given the (terminal) wind velocity, v_w.

Since the massive stars sampled in our simulations are typically hotter than the bi-stability-jump (with effective temperatures T_eff > 30 000 K), we used the fit formula of Eq. (3) in Vink & Sander (2021) to calculate v_w given the stellar luminosity, L, effective temperature, T_eff, and metallicity, Z. The wind velocity was then combined with L, T_eff, Z, the stellar mass, M, and escape velocity, $v_{\mathrm{esc}}=\sqrt{2GM/R_{\star }}$ (given the stellar radius R_⋆), to derive the mass-loss rate, ${\dot{m}}_{\mathrm{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_eff, and R_⋆ from T20 (EXT) to calculate v_w and ${\dot{m}}_{\mathrm{w}}$ for Z = 10⁻⁶, 10⁻⁵, 10⁻⁴, 0.01, and 1 Z_⊙. Here, we ignored the weak main-sequence stellar winds from extremely metal-poor stars with $Z<{10}^{-6}{\mathrm{Z}}_{\odot }$ 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, ${\dot{m}}_{\mathrm{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⁵¹ erg of thermal energy into their neighboring gas cells at the end of their main-sequence 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).

Since 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.

Kim & Ostriker (2015) show that the resolution is sufficiently 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

$$\begin{array}{c}\hfill M_{\mathrm{shell}}=1680{\mathrm{M}}_{\odot }{\left(\frac{E}{{10}^{51}\mathrm{erg}}\right)}^{0.87}{\left(\frac{n}{{\mathrm{cm}}^{-3}}\right)}^{-0.26}.\end{array}$$(12)

Following Equation (12), we can estimate the maximum density where the ST phase can be safely resolved, and it is 8 × 10⁶ cm⁻³ for our fiducial 1 M_⊙ resolution and 10³ cm⁻³ 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

$$\begin{array}{c}\hfill {\dot{N}}_{\mathrm{Ia}}(t)=N_{\mathrm{Ia}}{\displaystyle {\int }_0^t}{\dot{M}}_{\star }(t\prime )\mathrm{\Psi }(t-t\prime )\mathrm{d}t\prime ,\end{array}$$(13)

where N_Ia is an observable, ${\dot{M}}_{\star }(t)$ is the star formation rate (SFR), and Ψ(t) is the DTD. For each star particle in our simulations, ${\dot{M}}_{\star }(t)=m_{\star }\delta (t-t_0)$, where δ(t) is the Dirac Delta function, and t₀ is the formation time of the star particle. For simplicity, we denote t₀ = 0. This gives a simplified SN Ia rate for each star particle: ${\dot{N}}_{\mathrm{Ia}}(t)=m_{\star }{N}_{\mathrm{Ia}}\mathrm{\Psi }(t)$.

In 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, 𝒫_Ia, of being an SN Ia progenitor,

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{Ia}}=m_{\star }{N}_{\mathrm{Ia}}.\end{array}$$(14)

If 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⁻¹(y), where

$$\begin{array}{c}\hfill F(t)={\displaystyle {\int }_0^t}\mathrm{\Psi }(t\prime )\mathrm{d}t\prime .\end{array}$$(15)

We adopted a power-law DTD model that consists of a power law in time,

$$\begin{array}{c}\hfill \mathrm{\Psi }(t)=\{\begin{array}{ccc}& 0,\hfill & \hfill t<{\tau }_{\star }(8{\mathrm{M}}_{\odot },Z)\\ \hfill & A{t}^{-s},\hfill & \hfill t>{\tau }_{\star }(8{\mathrm{M}}_{\odot },Z),\end{array}\end{array}$$(16)

where A is a normalization factor such that the integration of Ψ(t) over a Hubble time, t_H, equals 1 (i.e., ∫₀^t_HΨ(t′)dt′ = 1) and τ_⋆(8 M_⊙, 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⁻³ SN ${\mathrm{M}}_{\odot }^{-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. Isolated dwarf galaxy simulations

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 MAKEDISK 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_{\mathrm{vir}}=2\times {10}^{10}{\mathrm{M}}_{\odot }$, 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⁷ M_⊙ and a scale length of 0.73 kpc and an initial gas disk with a mass of 4 × 10⁷ 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.

The gravitational softening lengths for dark matter (DM) and star were determined following Hopkins et al. (2018) by ${ϵ}_{\mathrm{DM}}=1.7\mathrm{pc}{m}_{\mathrm{DM}}^{1/3}{({\rho }_{\mathrm{DM}}/{10}^9{\mathrm{M}}_{\odot }{\mathrm{kpc}}^{-3})}^{-1/3}$ and ${ϵ}_{\star }=0.34\mathrm{pc}{m}_{\star }^{1/3}{(\langle n_{\mathrm{SF}}\rangle /{10}^3{\mathrm{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 $\stackrel{\sim }{c}=1000$ km s⁻¹ to capture the propagation of the ionization front in a typical ISM environment with a reasonable time lapse (Rosdahl et al. 2013).

To 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 Faucher-Giguè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⁴ cm⁻³ for 1 M_⊙ resolution and n_th = 10³ cm⁻³ 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.

In 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_⊙ is our fiducial run, named as “0.02 Z_⊙/full”. For comparison, we simulated the same galaxy without local radiative feedback (but with UVB), referred to as “0.02 Z_⊙/noRT”. We also simulated a galaxy initialized with 0.2 Z_⊙ with the “full” physics, and it is referred to as the “0.2 Z_⊙/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”.

We 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 different random seeds. We additionally ran the 0.02 Z_⊙/low simulation using another random seed, which demonstrates minor variations in the global SFR and ISM structures. Since different 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.

For each 1 M_⊙-resolution simulation, ∼1.3 × 10⁶ 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⁴ core-hours per 1 Gyr of simulation.

We 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 differences are minor.

We 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.

Fig. 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.02 Z⊙/noRT simulations show the most significant difference in cumulative stellar mass between the low- and high-resolution simulations; nonetheless, the difference remains within a factor of two.

The SFRs show an evident dependence on metallicity. The median SFR of 0.02 Z_⊙/full is only ∼6% of that in 0.2 Z_⊙/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⁴ 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 difficult for the ISM to revert to its cold phase once heated to its warm-diffuse phase due to ineffective cooling mechanisms and suppression of thermal instability (e.g. Schaye 2004; Walch et al. 2011; Bialy & Sternberg 2019).

Turning off radiative feedback has a larger impact on the SFR. The median SFR of 0.02 Z_⊙/noRT is 92 times higher than that of 0.02 Z_⊙/full and even five times higher than 0.2 Z_⊙/full. We will see in Section 3.4 that this difference is due to the fact that 0.02 Z_⊙/noRT contains the highest amount of cold and molecular gas, as there is no ISRF present to heat the diffuse gas and keep it warm and atomic.

We note that the SFRs exhibit good convergence with resolution. In all models, the cumulative stellar masses of the low-resolution simulations are slightly lower, but the difference in stellar mass formed between the low- and high-resolution simulations is within a factor of two. The 0.02 Z_⊙/noRT simulations show the most significant difference 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 differences 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).

To 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 high-resolution galaxies for the 0.02 Z_⊙/full (red), 0.2 Z_⊙/full (green), 0.02 Z_⊙/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.

Fig. 3. Star formation relations in the simulated dwarf galaxies. Left: Kennicutt-Schmidt relation for the simulated high-resolution galaxies for the 0.02 Z⊙/full (red), 0.2 Z⊙/full (green), and 0.02 Z⊙/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).

On 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.3 − 1.8 dex for the 0.15 < log(Σ_gas) < 0.75 range. The medians of 0.2 Z_⊙/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.02 Z_⊙/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 inefficient cooling. The 0.02 Z_⊙/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 efficient star formation process in the absence of momentum injection and photoheating from radiative feedback.

On the Σ_⋆^0.5Σ_gas − Σ_SFR plane (Fig. 3 middle panel), 0.2 Z_⊙/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.

The 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.2 Z_⊙/full galaxy, both the KS and extended KS relations are consistent with the observed empirical relations, while the results of 0.02 Z_⊙/full suggest a longer dynamical time-scale 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.02 Z_⊙/full, 0.2 Z_⊙/full, and 0.02 Z_⊙/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).

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 − Zinit; bottom row) of the simulated galaxies at 900 Myr. The three columns correspond to the different runs, as labeled above the panels. From left to right, they are 0.02 Z⊙/full, 0.2 Z⊙/full, and 0.02 Z⊙/noRT, respectively. The overset points represent young stars (< 100 Myr) with ages color-coded from violet (youngest) to red (oldest).

In 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.

The 0.02 Z_⊙/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.2 Z_⊙/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.02 Z_⊙/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 considerable 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 diffuse gas requires FUV radiation. Furthermore, the absence of radiation prevents efficient evacuation of the dense gas surrounding massive stars, which suppresses the efficiency of SN feedback. Consequently, the cold molecular gas can persist instead of being destroyed by FUV radiation that permeates the entire disc.

The 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.

Figure 5 presents the same projections and slices but viewed from an edge-on perspective. In the 0.02 Z_⊙/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.2 Z_⊙/full and 0.02 Z_⊙/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 CGM of the 0.02 Z_⊙/noRT galaxy is hotter and more enriched than that of the 0.2 Z_⊙/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.

Fig. 5. Similar to Fig. 4, but viewed from an edge-on perspective.

3.4. Interstellar medium properties and structures

3.4.1. Mass and volume fractions of different interstellar medium phases

Figure 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, 100 K < T < 8000 K), the warm ionized medium (WIM, 8000 K < T < 10 000 K), and hot ionized medium (HIM, T > 10 000 K) phases.

Fig. 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. Figure 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 efficient gas cooling and collapse (induced by higher cooling rate and lower heating rate in 0.2 Z⊙/full and 0.02 Z⊙/noRT, respectively), and the resultant higher SFRs.

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⁻³ of the ISM mass in the 0.02 Z_⊙/full galaxy, with its mass fraction increasing by a factor of around 10 in the 0.2 Z_⊙/full galaxy due to more effective metal cooling. Notably, in the 0.02 Z_⊙/noRT galaxy, the CNM mass fraction is about 40 times higher than in the 0.02 Z_⊙/full galaxy in the absence of radiative feedback. As a result, 0.2 Z_⊙/full and 0.02 Z_⊙/noRT exhibit higher SFRs due to the availability of more cold, dense gas for star formation. The higher SN rates in these two galaxies produce more hot gas as well. In 0.02 Z_⊙/full, the volume fraction of HIM is only ∼10⁻³ due to the low SN rate, while it constitutes a significant volume fraction (∼0.1) in 0.2 Z_⊙/full and 0.02 Z_⊙/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 differences with resolution, possibly due to the slightly lower SFR and potentially unresolved cooling at low 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 there 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.02 Z_⊙/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⁴ K isothermal line in 0.02 Z_⊙/full and 0.2 Z_⊙/full galaxies.

Fig. 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 G0 = 0.01, 0.1, and 0.00324 for 0.02 Z⊙/full, 0.2 Z⊙/full, and 0.02 Z⊙/noRT, respectively.

In 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 galaxies. 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⁴ K) and diffuse (∼1 cm⁻³). The distributions of warm gas in 0.02 Z_⊙/full and 0.2 Z_⊙/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⁻³, the gas in 0.2 Z_⊙/full is more susceptible to thermal instability, leading to easier cooling and condensation into cold, dense gas. Thus, 0.2 Z_⊙/full displays a distinct two-phase behavior with an additional peak of cold gas at temperatures ≲10² K. As a result, the 0.2 Z_⊙/full galaxy contains significantly more cold, dense gas than 0.02 Z_⊙/full. Additionally, the 0.2 Z_⊙/full also contains more hot, low-density gas because of its higher SN rate.

Fig. 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, diffuse 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 different 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 different metallicities and feedback channels are all dominated by the warm (∼104 K) diffuse (∼1 cm−3) H I gas. Radiative feedback leads to significantly less buildup of H2 in metal-poor environments.

Notably, the ISM phase structure of the 0.02 Z_⊙/noRT galaxy differs significantly from that of 0.02 Z_⊙/full, even for the warm diffuse 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 diffuse 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.02 Z_⊙/full, despite its lower metallicity and higher SFR. Furthermore, 0.02 Z_⊙/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).

Now, 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 H₂). In the second to fourth columns of Fig. 8, we present the mass distribution of gas for different galaxies weighted by the H II (purple), H I (orange), and H₂ (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⁻³ to at least 10³ cm⁻³, covering a broad range of temperatures of T ∼ 10 − 10⁴ 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⁻² cm⁻³), and by local photoionization in the vicinity of young massive stars. H II gas dominates the warm diffuse ISM and the hot low-density gas that corresponds to SN-driven bubbles, and it has an extended distribution towards the dense (∼10 cm³) cold (∼100 K) gas. In the 0.02 Z_⊙/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.2 Z_⊙/full, there is a departure from symmetry with a high-density tail comprising the photoionized gas within H II regions.

The H₂ 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, H₂ can only be destroyed slowly by collisional processes and CR ionization. In the 0.02 Z_⊙/full galaxy, H₂ remains a minority component even in dense gas with > 10³ cm³ due to its low metallicity. Half of the total H₂ mass exists in n_H < 1 cm⁻³ diffuse gas where the H₂ fraction is low. The situation is similar in the 0.2 Z_⊙/full galaxy. Although there is more cold, dense gas where H₂ can form, the H I–H₂ transition occurs at high densities of > 10³ cm⁻³ with a significant fraction of the H₂ mass existing in H I-dominated gas. The 0.02 Z_⊙/noRT galaxy exhibits a broader distribution of H₂ across large temperature and density ranges. H₂ even exists in warm-hot gas (T ∼ 10⁴ − 10⁶ K) with densities n_H ∼ 10⁻⁴ − 10⁻² cm⁻³, because the gas-phase formation of H₂ can be effective in the high-temperature diffuse 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 within the ISM regions of the two galaxies with radiative feedback, 0.02 Z_⊙/full and 0.2 Z_⊙/full. The ISFR energy density is described by the conventional parameter G₀ = u_{5.8 − 13.6 eV}/5.29 × 10⁻¹⁴ erg cm⁻³ (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.2 Z_⊙/full galaxy, the stacked ISRF peaks at around ∼0.1 G₀ in the central area and decreases exponentially to the background level at R ≃ 2 kpc. The ISRF in 0.02 Z_⊙/full is roughly one-tenth of that in the 0.2 Z_⊙/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.2 Z_⊙/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 different snapshots.

Fig. 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.02 Z⊙/full and 0.2 Z⊙/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 G0, UVB = 0.00 324.

The stacked ISRF in the low-resolution galaxies is also plotted with dotted curves, which show minor differences from the high-resolution results for the 0.2 Z_⊙ galaxies and the inner part of the 0.02 Z_⊙. However, the stacked ISRF in the 0.02 Z_⊙/full/low galaxy appears to be two to four times lower than that of 0.02 Z_⊙/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 efficiency 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.

Fig. 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.

In the case of the two full feedback galaxies, most of the SNe exploded in n_H < 1 cm⁻³ diffuse gas because of the preprocessing effect of radiation. The explosion density distributions of these two runs are similar at the high-density end, but for 0.2 Z_⊙/full, there are more SNe exploding in n_H < 0.01 cm⁻³ 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.2 Z_⊙/full, there are more clustered SNe because of its higher SFR, while most star clusters in 0.02 Z_⊙/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.2 Z_⊙/full.

In the absence of radiative feedback, the explosion density distribution of 0.02 Z_⊙/noRT is significantly broader, reaching up to 10⁴ cm⁻³. Although the main sequence winds can still provide early feedback, their ability to expel cold dense gas is limited by their low luminosities and efficient 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 inefficient in providing feedback and driving outflow (see Section 3.6 for the analysis of outflow properties).

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).

In vertical dynamical equilibrium of the galactic disk, the weight of the ISM per unit area, 𝒲, 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 𝒲, 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.

In 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_s²/(5/3), and the turbulent pressure P_turb = ρσ_z² are plotted separately. In the two full feedback galaxies, 0.02 Z_⊙/full and 0.2 Z_⊙/full, the ISM pressure is dominated by the thermal pressure which spans a similar range of 10^2.5 K cm⁻³ < P/k_B < 10^3.5 K cm⁻³. The turbulent pressure in these galaxies shows a wider variation, with P_turb being overall lower in 0.02 Z_⊙/full compared to 0.2 Z_⊙/full. The similarity in thermal pressure is due to the prevalence of the warm diffuse ISM, maintaining nearly identical thermal conditions when P/k_B < 10⁴ K cm⁻³. The lower turbulent pressure in 0.02 Z_⊙/full results from the reduced injection of turbulence attributed to its lower SFR. On the other hand, 0.02 Z_⊙/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.02 Z_⊙/full, attributed to its enhanced SFR.

Fig. 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.02 Z⊙/noRT.

Now 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 compared to 0.02 Z_⊙/noRT. The 0.02 Z_⊙/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.02 Z_⊙/full galaxy, showing that their fitting functions might still be applicable at metallicity lower than 0.1 Z_⊙. The difference 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 difference 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.

The feedback yields (P/Σ_SFR) of the two full feedback galaxies also have differences. The feedback yields for the total pressure are comparable when Σ_SFR is high ($\gtrsim {10}^{-4.5}{\mathrm{M}}_{\odot }$ yr⁻¹ kpc⁻²), but the 0.2 Z_⊙/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.02 Z_⊙/full, feedback is expected to yield higher thermal pressure compared to 0.2 Z_⊙/full, especially at lower Σ_SFR. This phenomenon can be attributed to inefficient 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.02 Z_⊙/full will yield a lower turbulent pressure than in 0.2 Z_⊙/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 low-metallicity environments results in a less efficient 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.

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_form–t_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 effects 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.2 Z_⊙/full galaxies, the CIMF follows a power-law distribution of $\mathrm{d}N/\mathrm{d}{M}_{\mathrm{cl}}\propto M_{\mathrm{cl}}^{-\alpha }$ with a slope α = 2.07, which is in agreement with the observed −2 slope (e.g. Portegies 2010; Krumholz et al. 2019). The CIMF of the 0.02 Z_⊙/full galaxies has a slightly steeper slope of −2.40, although the uncertainty of the slope is large due to low number statistics. Turning off 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 massive clusters and shaping the CIMF steeper. This indicates that radiative feedback disperses clouds quickly and efficiently.

Fig. 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 off the further star formation.

The power-law slopes for the low-resolution simulations are −2.16 for 0.2 Z_⊙/full/low, −2.34 for 0.02 Z_⊙/full/low, and −1.80 for 0.02 Z_⊙/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% to 90% 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.02 Z_⊙/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.

Fig. 13. Age spreads of the clusters (t90 − 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.

To 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 cluster (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 “+τ_⋆(100 M_⊙)” 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.02 Z_⊙/full and 0.2 Z_⊙/full, the time differences between t₁₀₀ 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 ineffective. Therefore, the clusters grow over long periods and accumulate higher final masses.

Fig. 14. Time between the formation of the first and last member star of a cluster (t100 − 0) as the function of the time of the first massive star in the cluster (tfms). The black dashed curves present the line of t100 − 0 = tfms + {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 t100 and tfms are smaller than 1 Myr for most clusters in the full feedback simulations, much shorter than the time required for the first SN explosion.

In 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⁴ M_⊙) young star clusters are closely correlated with the timing of feedback emergence. With effective 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 affect star formation activities, the ISM structure, 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 ${\dot{M}}_{\mathrm{out}}$, metal mass ${\dot{M}}_{\mathrm{Z},\mathrm{out}}$, momentum ${\dot{p}}_{\mathrm{out}}$, energy ${\dot{E}}_{\mathrm{out}}$) as the sum of q_k over all gas particles within a slab of thickness ΔL centered at height |z|

$$\begin{array}{c}\hfill {\dot{q}}_{\mathrm{out}}=\frac{1}{\mathrm{\Delta }L}{\displaystyle \sum _k}{v}_k·n{q}_k,(v_k·n>0),\end{array}$$(17)

where 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.

We 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 energy-loading (η_E) factors as

$$\begin{array}{c}\hfill {\eta }_M=\frac{{\dot{M}}_{\mathrm{out}}}{{\dot{M}}_{\star }},{\eta }_Z=\frac{{\dot{M}}_{\mathrm{Z},\mathrm{out}}}{Z_{\mathrm{ISM}}{\dot{M}}_{\star }},{\eta }_p=\frac{{\dot{p}}_{\mathrm{out}}}{p_{\mathrm{ej}}{\dot{N}}_{\mathrm{SN}}},{\eta }_E=\frac{{\dot{E}}_{\mathrm{out}}}{E_{\mathrm{SN}}{\dot{N}}_{\mathrm{SN}}},\end{array}$$(18)

respectively, where the SFR ${\dot{M}}_{\star }$ is time-averaged over 10 Myr, Z_ISM is the metallicity of the ISM region, ${\dot{N}}_{\mathrm{SN}}={\dot{M}}_{\star }/(100{\mathrm{M}}_{\odot })$ is the SN rate assuming a canonical IMF, p_ej = 3 × 10⁴ M_⊙ km s⁻¹ and E_SN = 10⁵¹ erg are the momentum carried by ejecta and injected energy of a typical SN.

In 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.02 Z_⊙/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.02 Z_⊙/full ones, due to their bursty star formation. On the contrary, the 0.02 Z_⊙/noRT results present the smallest fluctuations due to its more continuous star formation.

Fig. 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).

In 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 difference) loading factors for mass, metal, and energy, despite differences in their metallicities and SFRs. The momentum-loading factor of 0.02 Z_⊙/full is slightly higher than that of 0.2 Z_⊙/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 different in these two galaxies. The mass- and metal-loading factors of 0.02 Z_⊙/full and 0.2 Z_⊙/full fluctuate rapidly between 10 and 500, with a median of η_M ≈ 46, η_Z ≈ 45 for 0.02 Z_⊙/full and η_M ≈ 22 and η_Z ≈ 25 for 0.2 Z_⊙/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 difference 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 PdV work, while the total energy dissipates during the outflow process.

Remarkably, the mass-loading factor of 0.02 Z_⊙/noRT is lower than that of the full feedback simulations by almost an order 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 effect is mitigated in our simulations due to the low gas-surface density ($3{\mathrm{M}}_{\odot }{\mathrm{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 massive 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.02 Z_⊙/noRT, only 6% 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 effective 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 sufficient large clusters also determines the effectiveness of SN explosion in driving galactic outflows in the noRT case.

The 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.

4. Discussion

4.1. Effects of metallicity and radiative feedback on simulated dwarf galaxies

In this work, we simulated two galaxies alongside our fiducial simulation, one with ten times higher metallicity (0.2 Z_⊙/full) and another with radiative feedback turned off (0.02 Z_⊙/noRT), to study the effects of metallicity and radiative feedback.

The ten-times-higher metallicity results in a roughly ten-times-higher metal cooling rate in T > 10⁵ K hot gas and T < 10⁴ K warm 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.02 Z_⊙/full and 0.2 Z_⊙/full are consistent with observed dwarf galaxies (e.g. Roychowdhury et al. 2015) and exhibit superlinear relationships with similar slopes. The results of 0.2 Z_⊙/full roughly follow the empirical relations given by Elmegreen (2015) and Shi et al. (2018), while those of 0.02 Z_⊙/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 difference 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.2 Z_⊙/full because this galaxy can form larger molecular clouds that require more energetic feedback to disperse.

A 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 high-redshift 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.

On 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 effects 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⁻³ low-density gas. This effect increases the occurrence of efficient 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 high-resolution multiphase ISM simulations. The maximum gas density of the simulation is constrained by the Jeans criterion (m_gas > f_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_gas². In the absence of effective 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⁶ cm⁻³, the cool radius and cooling time of its ST phase are only 0.07 pc and 22 yr (Equations (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.

To 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 SN-only simulation is denoted as 0.02 Z_⊙/noESF, in comparison to 0.02 Z_⊙/noRT and 0.02 Z_⊙/full. Stars with extremely high metallicity (log(Z/Z_⊙) > 1) can be formed from the unmixed ejecta in 0.02 Z_⊙/noESF due to the absence of early feedback. Such stars have never been observed so far. In contrast, the stellar metallicity distribution in 0.02 Z_⊙/full exhibits a spread of ∼0.2 dex. While 0.02 Z_⊙/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.

Fig. 16. Stellar metallicity distribution of the simulation without any early feedback mechanism (black curve, 0.02 Z⊙/noESF) in comparison with 0.02 Z⊙/noRT and 0.02 Z⊙/full. Stars with extremely high metallicity (log(Z/Z⊙) > 1) can be formed from the unmixed ejecta in 0.02 Z⊙/noESF due to the absence of early feedback.

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.2 Z_⊙/full and 0.02 Z_⊙/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 self-regulating state. It is unknown if the SFR could further increase when more gas is available, leading to better agreement with our results.

The 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.2 Z_⊙/full galaxies, as reported by H17. In contrast, the volume fraction of HIM is only ∼10⁻³ in 0.02 Z_⊙/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 H₂ mass in warm diffuse gas with a low H₂ fraction in low-metallicity environments.

The ISRF of our 0.2 Z_⊙/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.

Similar 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.

The mass-loading factor of our 0.02 Z_⊙/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⁻³ so that they are more efficient to drive galactic outflow. This pre-processing effect of radiative feedback on the SN explosion site has been noticed in simulations of various environments (e.g., H17; Kannan et al. 2020b; 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 effects, 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.

As 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 difference 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 different 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 insufficient 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 efficiency, can vary significantly with the underlying sub-grid model (Hislop et al. 2022). Moreover, the IMF can be more top-heavy than the canonical one in metal-poor dwarf galaxies (e.g. Marks et al. 2012; Gennaro et al. 2018), thereby increasing the efficiency 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 efficiency (e.g., Hirai et al. 2021; L23).

Recent 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ć et al. (2021) for future study. We also do not model post-main-sequence 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 effects on the global properties 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 efficient 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 different 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.

We assume that all 8 − 100 M_⊙ massive stars have a constant SN energy of 10⁵¹ 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.

The metals diffuse passively following mesh advection, and we do not employ any sub-grid diffusion model for the unresolved eddies. In this case, the numerical diffusivity 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 diffusion 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 diffusion in detail will be a significant focus of the RIGEL project.

Dust formation and destruction in the gas are not explicitly accounted for in our current model. We assumed a constant dust-to-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.

Lastly, 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 different 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, affecting 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 low-mass reionization survivors (M_halo ≳ 10⁹ 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 star-by-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 efficient 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:

1. The RIGEL model self-consistently produces a multiphase ISM (Figs. 4 and 5) in which the star formation is regulated to a rate consistent with the observed dwarf galaxies (see Section 3.2).

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 Figs. 2 and 3) by removing the cold, dense gas fuel available for star formation (Fig. 6). Turning off radiative feedback boosts the SFR to an unreasonably high value by about two orders of magnitude.

3. The ISMs in three galaxies with different metallicities and feedback channels are all dominated by warm (100 K < T < 10⁵ K) H I gas in both mass and volume. Specifically, H I gas dominates at densities from n_H ∼ 10^−2.5 cm⁻³ to at least 10³ cm⁻³, covering a broad range of temperatures of T ∼ 10 − 10⁴ K. H II gas dominates the hot and photoionized gas. H₂ is always a minority phase in the 0.02 Z_⊙/full galaxy, even at a high density of n_H > 10³ cm⁻³ (see Fig. 8). With a higher metallicity, 0.2 Z_⊙/full is still dominated by H I below 10³ cm⁻³. Turning off radiative feedback leads to diffuse H₂ formed in warm-hot gas with temperatures and densities extended to 10⁶ K and 10⁻⁴ cm⁻³ in the absence of the LW radiation.

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.2 Z_⊙/full is roughly ten times higher than 0.02 Z_⊙/full because of its higher SFR.

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 efficient 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 off all the early feedback results in a rapid cooling of (unmixed) SN ejecta within dense gas, leading to the formation of excessively metal-rich stars.

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 effect on the ISM (see Section 3.4.5).

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 efficient 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).

8. The mass-loading factors of the galactic outflows in our full feedback galaxies vary between 10 and 500 (Fig. 15). Turning off the radiative feedback leads to reduced mass- and momentum-loading factors, while having minor effects 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.

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.

Acknowledgments

We thank 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 NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), ASTROPY (Astropy Collaboration 2013, 2018), and MATPLOTLIB (Hunter 2007) to analyze and visualize the simulation data.

References

Appendix A: Treatments of cooling channels

In this appendix, we provide a detailed description of the various cooling channels in our model.

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. (2020a). We adopt the rate coefficients accounting for cooling/heating of atomic hydrogen and helium from Cen (1992). We model the cooling due to molecular hydrogen by

$$\begin{array}{c}\hfill n^2{\mathrm{\Lambda }}_{{\mathrm{H}}_2}={\mathrm{\Lambda }}_{\text{H}\text{I}{\mathrm{H}}_2,n\to 0}(T)n_{\text{H}\text{I}}{n}_{{\mathrm{H}}_2}+{\mathrm{\Lambda }}_{{\mathrm{H}}_2{\mathrm{H}}_2,n\to 0}(T)n_{{\mathrm{H}}_2}{n}_{{\mathrm{H}}_2},\end{array}$$(A.1)

where ${\mathrm{\Lambda }}_{\mathrm{H}\mathrm{I}{\mathrm{H}}_2,n\to 0}$ and Λ_{H2H2, n → 0} are the H I–H₂ and H₂–H₂ collisional cooling coefficients in the low-density limit (n → 0) from Hollenbach & McKee (1979). We temporarily ignore UV pumping heating and H₂ formation heating as they are subdominant in our current framework.

A.2. Metal collisional excitation lines in hot gas

The metal-line cooling in ≳10⁵ 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 self-shielding 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).

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:

$$\begin{array}{cc}\hfill n^2{\mathrm{\Lambda }}_{\mathrm{neb}}& =3.86\times {10}^{-23}{\mathrm{erg}\mathrm{cm}}^{-3}{\mathrm{s}}^{-1}\hfill \\ \hfill & \times Z_{\mathrm{g}}\prime n_e{n}_{\text{H}\text{II}}\frac{e^{-3.86/T_4}{f}_{\mathrm{neb}}(T)}{T_4^{1/2}(1+0.12n_{e,2}^{0.38-0.12ln{T}_4})},\hfill \end{array}$$(A.2)

where Z_g′ is the gas metallicity normalized to the solar value, T₄ = T/10⁴ K, n_e, 2 = n_e/10² cm⁻³, and ${log}_{10}{f}_{\mathrm{neb}}(T)={\displaystyle \sum _{i=0}^6}{a}_i{(ln{T}_4)}^i$ with {a₀, a₁, a₂, a₃, a₄, a₅}={0.692, −0.586, 0.816, −0.505, 0.118, 7.66 × 10⁻³, −5.08 × 10⁻³}. 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).

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 other species (e.g. Wolfire et al. 1995, 2003). As noticed by Kim et al. (2023a), the low-temperature (< 10⁴ 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 coefficients given by Gong et al. (2017). For example, the cooling rate for the C II 158 μm two-level system is

$$\begin{array}{c}\hfill n^2{\mathrm{\Lambda }}_{\text{C}\text{II}}=n_{\text{C}\text{II}}{A}_{\mathrm{ul}}{E}_{\mathrm{ul}}\frac{n_j{k}_{\mathrm{lu}}}{n_j{k}_{\mathrm{lu}}+n_j{k}_{\mathrm{ul}}+A_{\mathrm{ul}}},\end{array}$$(A.3)

where A_ul is the Einstein A coefficient, E_ul is the energy difference, k_ul is the collisional rate coefficien and k_lu = g_u/g_lk_ulexp( − E_ul/k_BT), and n_j is the number density of collisional partners including H₂, 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).

Fig. A.1. Cooling functions in predominantly neutral gas with nH = 10 cm−3, xe = 2.7 × 10−3, G0 = 1 for different 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.

In Fig. A.1, we compare the cooling functions in predominantly neutral gas with n_H = 10 cm⁻³, x_e = 2.7 × 10⁻³, G₀ = 1 obtained 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 difference. 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 43n_H) in their fitting function.

Compared to FIRE-3, RIGEL demonstrates overall higher cooling rates. This variation occurs due to two main factors. First, RIGEL incorporates the cooling effects 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 coefficients 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.

A.5. CO cooling in cold gas

In dense (n ≳ 10³ cm⁻³) 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 & Jaffa (2018),

$$\begin{array}{c}\hfill n^2{\mathrm{\Lambda }}_{\mathrm{CO}}=n_{\mathrm{CO}}{x}_{{\mathrm{H}}_2}{({\mathrm{\Lambda }}_{\mathrm{CO},\mathrm{LO}}^{-1/\beta }+{\mathrm{\Lambda }}_{\mathrm{CO},\mathrm{HI}}^{-1/\beta })}^{-\beta },\end{array}$$(A.4)

where Λ_CO, LO, Λ_CO, HI, and β are given by

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_{\mathrm{CO},\mathrm{LO}}& =2.16\times {10}^{-27}{\mathrm{erg}\mathrm{s}}^{-1}\left(\frac{n_{\mathrm{H}2}}{{\mathrm{cm}}^{-3}}\right){\left(\frac{T}{K}\right)}^{3/2},\hfill \\ \hfill {\mathrm{\Lambda }}_{\mathrm{CO},\mathrm{HI}}& =2.21\times {10}^{-28}{\mathrm{erg}\mathrm{s}}^{-1}\hfill \\ \hfill & \times {\left[\frac{(x_{\mathrm{CO}}/x_{{\mathrm{H}}_2}){\mathrm{km}\mathrm{s}}^{-1}{\mathrm{pc}}^{-1}}{|\mathrm{\nabla }·\mathit{v}|}\right]}^{-1}{\left(\frac{n_{{\mathrm{H}}_2}}{{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\frac{T}{K}\right)}^4,\hfill \\ \hfill \beta & =1.23{\left(\frac{n_{{\mathrm{H}}_2}}{{\mathrm{cm}}^{-3}}\right)}^{0.0533}{\left(\frac{T}{K}\right)}^{0.164}.\hfill \end{array}$$(A.5)

which reproduces the detailed numerical results of Goldsmith & Langer (1978) for T ∈ [10,100] K, n_H2 ∈ [3×10²,3×10⁵] cm⁻³, and |∇⋅v|[x_CO/x_H2/(3 × 10⁻⁴)]⁻¹ ∈ [10,10⁶] km s⁻¹ pc⁻¹.

A.6. Dust cooling in cold gas

In even denser molecular gas (n ≳ 10⁶ cm⁻³) 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. (2020a):

$$\begin{array}{cc}\hfill n^2{\mathrm{\Lambda }}_{\mathrm{gd}}& =1.356\times {10}^{-33}{\mathrm{erg}\mathrm{cm}}^{-3}{\mathrm{s}}^{-1}\hfill \\ \hfill & \times \left(\frac{Z_{\mathrm{d}}}{0.01}\right)\left(\frac{0.1\mu m}{a}\right)\sqrt{\frac{T}{\mathrm{K}}}\left(\frac{T-T_{\mathrm{d}}}{\mathrm{K}}\right){\left(\frac{n_{\mathrm{H}}}{{\mathrm{cm}}^{-3}}\right)}^2.\hfill \end{array}$$(A.6)

where 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 instantaneous 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

$$\begin{array}{c}\hfill n^2{\mathrm{\Gamma }}_{\mathrm{dr}}={\kappa }_{\mathrm{P}}\rho (4\sigma T_{\mathrm{d}}^4-c{E}_{\mathrm{IR}}),\end{array}$$(A.7)

where κ_P is the Planck mean opacity, σ is the Stefan-Boltzmann constant, and E_IR is the energy density of IR photons.

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. (2020a), we adopt the photoelectric heating rate given by Wolfire et al. (2003)

$$\begin{array}{c}\hfill n^2{\mathrm{\Gamma }}_{\mathrm{PE}}=1.3\times {10}^{-24}{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{cm}}^{-3}\left(\frac{Z_{\mathrm{d}}}{0.01}\right){ϵ}_{\mathrm{ff}}{G}_0{n}_{\mathrm{H}},\end{array}$$(A.8)

where G₀ is the FUV radiation energy density relative to Habing (1968)’s estimate (u_{5.8 − 13.6eV} = 5.29 × 10⁻¹⁴ erg cm⁻³) and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {ϵ}_{\mathrm{ff}}& =\frac{0.0487}{1+0.004{[G_0\sqrt{T}/(0.5n_e)]}^{0.73}}\hfill \\ \hfill & +\frac{0.0365{(T/{10}^4)}^{0.7}}{1+2\times {10}^{-4}\times G_0\sqrt{T}/(0.5n_e)}.\hfill \end{array}\end{array}$$(A.9)

Here, the unit of temperature T is [K], and the unit of electron number density n_e is [cm⁻³].

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

$$\begin{array}{cc}\hfill {\zeta }_{\mathrm{cr},{\mathrm{H}}_2}& =3.5\times {10}^{-16}{\mathrm{s}}^{-1}\hfill \\ \hfill {\zeta }_{\mathrm{cr},\text{H}\text{I}}& =1.78\times {10}^{-16}{\mathrm{s}}^{-1}\hfill \\ \hfill {\zeta }_{\mathrm{cr},\text{He}\text{I}}& =1.1{\zeta }_{\mathrm{cr},}\text{H}\text{I}\hfill \\ \hfill {\zeta }_{\mathrm{cr},\text{C}\text{I}}& =3.85{\zeta }_{\mathrm{cr},}\text{H}\text{I}\hfill \\ \hfill n^2{\mathrm{\Gamma }}_{\mathrm{cr},j}& =1.6022\times {10}^{-11}\mathrm{erg}·{\zeta }_{\mathrm{cr},j}{n}_j,\hfill \end{array}$$(A.10)

where $j\in {\mathrm{H}}_2,\mathrm{H}\mathrm{I},\mathrm{He}\mathrm{I}$.

Appendix B: Tables of fitting parameters

In Table B.1, Table B.2, and Table B.3, we listed the fitting parameters x₀ and b_i(i = 0, 1, 2, 3, 4) the for main-sequence lifetime (in [yr]), UV photon production rates in different bands (in [s⁻¹]), and wind mass-loss rate (in $[{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}]$) and velocity (in [km s⁻¹]), defined in Eq. (8), respectively. The machine-readable tables with full precision are downloadable online⁵.

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).

As 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 low-density gas (Section 3.4.4), where the 10 M_⊙ resolution is sufficient to resolve all SNe. In the 0.02Z_⊙/noRT/low simulation, a small fraction of SNe is marginally resolved, while they should have minor effects on the results.

Fig. C.1. Cumulative distribution function of the ratio between the mass of the cell where the SNe injected explosion energy and the total swept-up 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 ST-resolving criterion of mgas, SN / Mshell < 1/27 by Kim & Ostriker (2015).

All Tables

All Figures

Fig. 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}{\mathrm{Z}}_{\odot }$) 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}{\mathrm{Z}}_{\odot }$) to high (1 Z⊙). We shifted the curves of the Opt., FUV, LW, and EUV1 bands by multiplying constants of 104, 103, 102, and 10, respectively, to avoid overlapping different bins. Right: Similar to left but for the wind mass-loss rate (${\dot{m}}_{\mathrm{w}}$) and wind velocity (vw). The color of the curves from light to dark indicates the metal abundance from low (${10}^{-6}{\mathrm{Z}}_{\odot }$) to high (1 Z⊙). The main-sequence stellar winds from extremely metal-poor stars with $Z<{10}^{-6}{\mathrm{Z}}_{\odot }$ are negligible and not considered in our model.

In the text

Fig. 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.02 Z⊙/noRT simulations show the most significant difference in cumulative stellar mass between the low- and high-resolution simulations; nonetheless, the difference remains within a factor of two.

In the text

Fig. 3. Star formation relations in the simulated dwarf galaxies. Left: Kennicutt-Schmidt relation for the simulated high-resolution galaxies for the 0.02 Z⊙/full (red), 0.2 Z⊙/full (green), and 0.02 Z⊙/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).

In the text

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 − Zinit; bottom row) of the simulated galaxies at 900 Myr. The three columns correspond to the different runs, as labeled above the panels. From left to right, they are 0.02 Z⊙/full, 0.2 Z⊙/full, and 0.02 Z⊙/noRT, respectively. The overset points represent young stars (< 100 Myr) with ages color-coded from violet (youngest) to red (oldest).

In the text

	Fig. 5. Similar to Fig. 4, but viewed from an edge-on perspective.
In the text

Fig. 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. Figure 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 efficient gas cooling and collapse (induced by higher cooling rate and lower heating rate in 0.2 Z⊙/full and 0.02 Z⊙/noRT, respectively), and the resultant higher SFRs.

In the text

	Fig. 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 G0 = 0.01, 0.1, and 0.00324 for 0.02 Z⊙/full, 0.2 Z⊙/full, and 0.02 Z⊙/noRT, respectively.
In the text

Fig. 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, diffuse 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 different 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 different metallicities and feedback channels are all dominated by the warm (∼104 K) diffuse (∼1 cm−3) H I gas. Radiative feedback leads to significantly less buildup of H2 in metal-poor environments.

In the text

Fig. 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.02 Z⊙/full and 0.2 Z⊙/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 G0, UVB = 0.00 324.

In the text

	Fig. 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.
In the text

Fig. 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.02 Z⊙/noRT.

In the text

Fig. 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 off the further star formation.

In the text

Fig. 13. Age spreads of the clusters (t90 − 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.

In the text

Fig. 14. Time between the formation of the first and last member star of a cluster (t100 − 0) as the function of the time of the first massive star in the cluster (tfms). The black dashed curves present the line of t100 − 0 = tfms + {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 t100 and tfms are smaller than 1 Myr for most clusters in the full feedback simulations, much shorter than the time required for the first SN explosion.

In the text

Fig. 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).

In the text

	Fig. 16. Stellar metallicity distribution of the simulation without any early feedback mechanism (black curve, 0.02 Z⊙/noESF) in comparison with 0.02 Z⊙/noRT and 0.02 Z⊙/full. Stars with extremely high metallicity (log(Z/Z⊙) > 1) can be formed from the unmixed ejecta in 0.02 Z⊙/noESF due to the absence of early feedback.
In the text