# Metallicity Dependence of Pressure-regulated Feedback-modulated Star Formation in the TIGRESS-NCR Simulation Suite

Chang-Goo Kim, Eve C. Ostriker, Jeong-Gyu Kim, Munan Gong, Greg L.

Bryan, Drummond B. Fielding, Sultan Hassan, Matthew Ho, Sarah M. R.

Jeffreson, Rachel S. Somerville, et al.

► To cite this version:

Chang-Goo Kim, Eve C. Ostriker, Jeong-Gyu Kim, Munan Gong, Greg L. Bryan, et al. Metallicity Dependence of Pressure-regulated Feedback-modulated Star Formation in the TIGRESS-NCR Simulation Suite. The Astrophysical Journal, 2024, 972, 10.3847/1538-4357/ad59ab. insu-04851533

HAL Id: insu-04851533

https://insu.hal.science/insu-04851533v1

Submitted on 20 Dec 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Metallicity Dependence of Pressure-regulated Feedback-modulated Star Formation in the TIGRESS-NCR Simulation Suite

Chang-Goo Kim<fnref n="1" />, Eve C. Ostriker<fnref n="1" />,<fnref n="2" />, Jeong-Gyu Kim<fnref n="3" />,<fnref n="4" />, Munan Gong<fnref n="5" />, Greg L. Bryan<fnref n="6" />, Drummond B. Fielding<fnref n="7" />,<fnref n="8" />, Sultan Hassan<fnref n="7" />,<fnref n="9" />,<fnref n="10" />,<fnref n="13" />, Matthew Ho<fnref n="11" />, Sarah M. R. Jeffreson<fnref n="12" />, Rachel S. Somerville<fnref n="7" />, and Ulrich P. Steinwandel<fnref n="7" />

Received 2024 April 3; revised 2024 May 28; accepted 2024 June 17; published 2024 August 26

## Abstract

We present a new suite of numerical simulations of the star-forming interstellar medium (ISM) in galactic disks using the TIGRESS-NCR framework. Distinctive aspects of our simulation suite are (1) sophisticated and comprehensive numerical treatments of essential physical processes including magnetohydrodynamics, self-gravity, and galactic differential rotation, as well as photochemistry, cooling, and heating coupled with direct ray-tracing UV radiation transfer and resolved supernova feedback and (2) wide parameter coverage including the variation in metallicity over $Z' \equiv Z/Z_\odot \sim 0.1-3$, gas surface density $\Sigma_{\text{gas}} \sim 5-150 M_\odot \text{pc}^{-2}$, and stellar surface density $\Sigma_{\text{star}} \sim 1-50 M_\odot \text{pc}^{-2}$. The range of emergent star formation rate surface density is $\Sigma_{\text{SFR}} \sim 10^{-4}-0.5 M_\odot \text{kpc}^{-2} \text{yr}^{-1}$, and ISM total midplane pressure is $P_{\text{tot}}/k_B = 10^3-10^6 \text{cm}^{-3} \text{K}$, with $P_{\text{tot}}$ equal to the ISM weight $W$. For given $\Sigma_{\text{gas}}$ and $\Sigma_{\text{star}}$, we find $\Sigma_{\text{SFR}} \propto Z'^{0.3}$. We provide an interpretation based on the pressure-regulated feedback-modulated (PRFM) star formation theory. The total midplane pressure consists of thermal, turbulent, and magnetic stresses. We characterize feedback modulation in terms of the yield $\Upsilon$, defined as the ratio of each stress to $\Sigma_{\text{SFR}}$. The thermal feedback yield varies sensitively with both weight and metallicity as $\Upsilon_{\text{th}} \propto W^{-0.46}Z'^{-0.53}$, while the combined turbulent and magnetic feedback yield shows weaker dependence $\Upsilon_{\text{turb+mag}} \propto W^{-0.22}Z'^{-0.18}$. The reduction in $\Sigma_{\text{SFR}}$ at low metallicity is due mainly to enhanced thermal feedback yield, resulting from reduced attenuation of UV radiation. With the metallicity-dependent calibrations we provide, PRFM theory can be used for a new subgrid star formation prescription in cosmological simulations where the ISM is unresolved.

*Unified Astronomy Thesaurus concepts: Interstellar medium (847); Star formation (1569); Magnetohydrodynamical simulations (1966); Stellar feedback (1602); Metallicity (1031); Galaxy formation (595); Radiative transfer simulations (1967)*

## 1. Introduction

Galactic star formation rates (SFRs) and the physical state of the interstellar medium (ISM) are observed to be tightly connected (e.g., Kennicutt & Evans 2012; Sun et al. 2020, 2023; Barrera-Ballesteros et al. 2021). This connection can be understood theoretically based on ISM dynamics and thermodynamics and the physics of stellar feedback. On the one hand, the loss of energy in the ISM occurs on relatively short timescales via radiative cooling and turbulence dissipation. Gas with locally reduced pressure support in turn collapses by gravity and forms stars. On the other hand, newborn stars return energy (sourced by nuclear fusion) that pervades the surrounding ISM, offsetting losses and recovering

the balance between pressure (provided by thermal, turbulent, and magnetic components) and gravity (e.g., Ostriker et al. 2010; Kim et al. 2011; Ostriker & Shetty 2011; Ostriker & Kim 2022). Stellar feedback is thus a key in controlling future star formation and maintaining the physical state of ISM disks. Because ISM evolution, star formation, and feedback are inherently cyclic, a realistic understanding of the galactic ecosystem necessitates a holistic approach to these tightly coupled physical processes.

Representing ISM physics and stellar feedback in numerical simulations requires treatments of the thermodynamic properties of gas in different phases and of localized injection of energy from stellar feedback. With varying degrees of accuracy, there exist many numerical frameworks that solve (magneto)hydrodynamics equations including losses—from shocks, turbulent cascades, and phase mixing followed by radiative cooling—and modeling gains—energy returns mainly from massive young stars—in the context of the galactic ecosystem. Such efforts can be categorized into three different types based on their outer scales: cosmological zoom-in

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA; cgkim@astro.princeton.edu</footnote>

<footnote n="2">²Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA</footnote>

<footnote n="3">³Division of Science, National Astronomical Observatory of Japan, Mitaka, Tokyo 181-0015, Japan</footnote>

<footnote n="4">⁴Korea Astronomy and Space Science Institute, Daejeon 34055, Republic of Korea</footnote>

<footnote n="5">⁵Max-Planck Institute for Extraterrestrial Physics, Garching near Munich, D-85748, Germany</footnote>

<footnote n="6">⁶Department of Astronomy, Columbia University, 550 West 120th Street, New York, NY 10027, USA</footnote>

<footnote n="7">⁷Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA</footnote>

<footnote n="8">⁸Department of Astronomy, Cornell University, Ithaca, NY 14853, USA</footnote>

<footnote n="9">⁹Center for Cosmology and Particle Physics, Department of Physics, New York University, 726 Broadway, New York, NY 10003, USA</footnote>

<footnote n="10">¹⁰Department of Physics & Astronomy, University of the Western Cape, Cape Town 7535, South Africa</footnote>

<footnote n="11">¹¹CNRS & Sorbonne Université, Institut d'Astrophysique de Paris (IAP), UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="12">¹²Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="13">¹³NASA Hubble fellow.</footnote>
simulations (e.g., Hopkins et al. 2014, 2018, 2023; Buck et al. 2020), isolated global galaxy simulations (e.g., Hu et al. 2017, 2023a, 2023b; Kannan et al. 2020b; Li et al. 2020; Smith et al. 2020, 2021; Jeffreson et al. 2021; Treß et al. 2021; Bieri et al. 2023; Steinwandel et al. 2023, 2024; Li et al. 2024), and vertically stratified local simulations of galactic disks (e.g., Gatto et al. 2017; Kim & Ostriker 2017; Kannan et al. 2020a; Kim et al. 2020b, 2023a; Hu et al. 2021, 2023b; Rathjen et al. 2021, 2023). Given the limited dynamic range of any numerical simulation, larger outer scales also imply coarser resolution of the ISM, which makes it difficult to follow multiphase physics explicitly. In particular, following the creation and evolution of the hot ISM is challenging for pseudo-Lagrangian approaches because it is so diffuse; at realistic hot densities $n_H < 10^{-2} \text{cm}^{-3}$, even a $(10 \text{ pc})^3$ volume contains $<0.3 \text{ M}_\odot$. Limited resolution also means that when gravitational collapse occurs, the mass involved may be more strongly clustered than is realistic, leading to an excessive spatio-temporal correlation of feedback. To date, larger-scale simulations of cosmological volumes typically have not attempted to explicitly resolve the multiphase ISM due to their insufficient resolution, instead adopting subgrid models such as those of Springel & Hernquist (2003) or Schaye & Dalla Vecchia (2008) to model star formation (e.g., Vogelsberger et al. 2020). However, efforts to move toward more explicit treatments of key physics in cosmological simulations are underway (e.g., Dubois et al. 2021; Feldmann et al. 2023).

Localized energy injection from stellar feedback results in expanding bubbles of different astronomical types, which have been studied using targeted numerical simulations. These include supernova remnants (SNRs) and superbubbles (Kim & Ostriker 2015a; Iffrig & Hennebelle 2015; Martizzi et al. 2015; Walch & Naab 2015; Kim et al. 2017a; Fielding et al. 2018; El-Badry et al. 2019; Gentry et al. 2019; Steinwandel et al. 2020), H II regions (Geen et al. 2016, 2020; Kim et al. 2018, 2021; Deng et al. 2024), and stellar wind blown bubbles (Haid et al. 2018; Lancaster et al. 2021a, 2021b; Geen et al. 2021); regions with strong feedback interact with each other and the surrounding inhomogeneous ISM. Expanding feedback-driven bubbles inject a significant amount of radial momentum in the ISM and cause phase transitions to hotter phases by photoionization and shocks in the gas (see, e.g., reviews of Ostriker & McKee 1988; Krumholz et al. 2019; Chevance et al. 2023). It is also important to note that the motions driven by feedback are coupled with other large-scale flows in galactic disks, e.g., shear and epicyclic motions induced by galactic differential rotation, and other externally driven gas flows, e.g., cosmic accretion.

Each of the different types of holistic numerical simulations has advantages and disadvantages. Cosmological zoom-in simulations can realistically capture the cosmic flows that build a given galaxy and interactions with other galaxies, while sacrificing accuracy in modeling the multiphase gas. In order to model feedback at low resolution in cosmological zoom simulations, the current best practice is to inject the terminal radial momentum that feedback bubbles should have (mainly due to supernovae, hereafter SNe) as calibrated from higher-resolution simulations (e.g., Hopkins et al. 2014; Kimm & Cen 2014; Oku et al. 2022). While this approach can drive turbulence in warm and cold ISM phases, it misses transitions to the hot phase. Doing so would require much higher resolution to follow the energy-conserving stage of SNRs and

shocks that create hot gas (e.g., Kim & Ostriker 2015a; Steinwandel et al. 2020). In cosmological zoom and global galaxy simulations, treatments of thermodynamic, chemical, and radiative processes in the ISM also involve many approximations. Local tall-box simulations as cited above are in the opposite limit, in which the ISM physics are treated with more accurate and explicit methods, and SN shock heating can be directly resolved.<fnref n="14" /> Some local (or semiglobal) models include galactic differential rotation using a shearing box approximation (Kim & Ostriker 2017; Colling et al. 2018) and galactic structures like spiral arms (Kim et al. 2020c) and bar driven inflows (Moon et al. 2021, 2023). However, the effects of global geometry and cosmic inflows cannot be captured directly. Isolated global galaxy simulations are in between; models of more massive galaxies (Milky Way–like) tend to have numerical approaches closer to cosmological zooms (e.g., Jeffreson et al. 2021, 2022; Lee et al. 2022, which treats feedback via prescribed momentum injection), while models of less massive galaxies (dwarf galaxies) include more explicit ISM physics (Kannan et al. 2020b; Katz et al. 2022; Katz 2022) and directly resolve SN feedback (Hu et al. 2023b; Steinwandel & Goldberg 2023; Steinwandel et al. 2023, 2024).

In a recent publication (Kim et al. 2023a), we presented the first results from simulations employing the new TIGRESS-NCR<fnref n="15" /> framework to study the star-forming ISM in conditions similar to those in the solar neighborhood and in local-universe galaxies such as those observed by PHANGS (Schinnerer et al. 2019). Here, we shall present results from a much larger simulation survey conducted using TIGRESS-NCR. The range of gas and stellar surface density in our new parameter survey is similar to that explored in Kim et al. (2020b) and Ostriker & Kim (2022) using the original TIGRESS framework (Kim & Ostriker 2017; we shall refer to this framework as “TIGRESS-classic” henceforth). The TIGRESS-classic framework solves ideal MHD equations within a local shearing box, utilizing uniformly high resolution (~2–8 pc), which enables gravitational collapse in dense regions—producing sink particles that act as feedback sources, while also being sufficient to resolve the energy-conserving stage of SNRs.

TIGRESS-NCR extends the TIGRESS-classic framework by including explicit UV radiation transfer using an adaptive ray-tracing (ART) method (Kim et al. 2017b) coupled with nonequilibrium photochemistry, cooling, and heating, as detailed in Kim et al. (2023b). This comprehensive “NCR” treatment of microphysics replaces the simplified treatment in TIGRESS-classic, which adopted approximate formulae for heating and cooling in warm-cold gas and was restricted to solar metallicity. Taking advantage of our expanded ISM modeling capabilities, the new suite of TIGRESS-NCR simulations covers from supersolar to low-metallicity regimes.

<footnote n="14">¹⁴ We note that there is an additional class of local simulations that do not follow self-consistent cycles of star formation and feedback (e.g., Joung & Mac Low 2006; Joung et al. 2009; Walch et al. 2015; Girichidis et al. 2016, 2018; Martizzi et al. 2016; Fielding et al. 2018; Tan & Fielding 2024). Rather, these models focus on the response of the ISM to prescribed injection of energy mimicking stellar feedback and use the local box setup in the interest of more control and higher resolution (see Schneider et al. 2020; Schneider & Mao 2024 for a similar example of global galaxy models). These numerical experiments have specific goals and merits, but work of this kind should not be confused with self-consistent simulations of the star-forming ISM like those presented in this paper, which aim to provide a holistic view of coregulation of the ISM and star formation.</footnote>

<footnote n="15">¹⁵ “TIGRESS” stands for “Three-phase ISM in Galaxies Resolving Evolution with Star formation and Stellar feedback,” and “NCR” stands for “None-equilibrium Cooling and Radiation.”</footnote>
Combining more than a decade variation in metallicity with a wide range of galactic conditions (gas surface density and gravitational potential), the simulation parameter study presented here allows us to broadly characterize scaling relations of galactic SFRs and turbulence. We note that effects of varying metallicity were previously studied in local-box star-forming ISM simulations by Hu et al. (2021). They explored the same metallicity range as we do, for solar neighborhood conditions, but the main scientific focus of that paper was on the metallicity dependence of the atomic-to-molecular transition and C⁺/C/CO distributions.

The exploration of metallicity dependence is critical as the thermal balance in the ISM is sensitive to the abundances of metals and dust.<fnref n="16" /> A variety of processes are responsible for radiative cooling and heating in different gas phases (Draine 2011; Ferland et al. 2017; Wolfire et al. 2022; Kim et al. 2023b). In the cold neutral ISM, fine-structure metal lines from C⁺ and O are the main cooling channels in atomic and diffuse molecular gas, with rotational transitions of CO becoming the dominant coolant in denser and more shielded gas; far-UV (FUV) producing the photoelectric (PE) effect in small grains and polycyclic aromatic hydrocarbons (PAHs) and cosmic-ray (CR) ionization are the main heating channels (Field et al. 1969; Watson 1972; Bakes & Tielens 1994; Wolfire et al. 1995, 2003; Weingartner & Draine 2001a; Bialy & Sternberg 2019). In the warm neutral ISM, cooling is dominated by the hydrogen Lyα line, while PE heating still dominates when grain abundances are sufficiently high; CR heating begins to dominate at low dust abundance. In the warm ionized ISM (Haffner et al. 2009), extreme-UV (EUV) photons ionize hydrogen, helium, and heavier atoms, making photo-ionization heating the dominant heating channel, while photoionized metals (e.g., O⁺, O⁺², N⁺) are the major coolants (we collectively refer to this as nebular line cooling; see Ferland et al. 2017 for comprehensive nebular modeling). Finally, hot gas ($T > 10^6$ K) is created by very high-speed shocks generated by SNe (Cox 1972; McKee & Ostriker 1977) and cools very inefficiently. Shocks at somewhat lower speed and mixing of the hot gas with denser material lead to temperatures $T \sim 10^{5-6}$ K, at which highly ionized metals become the dominant coolants (Sutherland & Dopita 1993; Gnat & Sternberg 2007; Gnat & Ferland 2012).

Metallicity affects ISM thermodynamics not only through direct agents of cooling (with lower metals reducing cooling) and heating (with lower dust reducing heating) but also through effects on UV radiative transfer (see Wolfire et al. 2022, for a review). Since grains absorb both FUV and EUV photons as they propagate through the ISM, lower abundances of dust reduce the attenuation of radiation and thereby tend to increase the gas heating for a given rate of radiation production. Because of the complex interplay among the different effects involved, a quantitative understanding of ISM thermodynamics at varying metallicity requires numerical simulations.

In the pressure-regulated, feedback-modulated (PRFM) theory of the star-forming ISM (see Ostriker & Kim 2022, and references therein), the ISM pressure varies directly with the SFR per unit area because feedback is responsible for heating and driving of turbulence. This relationship is

quantified in terms of the feedback yield $\Upsilon$ (see Equation (13) below for definition), which has previously been measured using TIGRESS-classic simulations in Ostriker & Kim (2022) and TIGRESS-NCR simulations in Kim et al. (2023a). Since midplane pressure is regulated to match the ISM weight under vertical equilibrium in disk galaxies, the feedback yield can be used to predict the mean SFR given large-scale galactic properties (primarily gas and stellar surface density); this prediction, as well as predictions for the relations between component pressures and star formation, has been validated in nearby galaxies (Herrera-Camus et al. 2017; Sun et al. 2020, 2023; Barrera-Ballesteros et al. 2021). A question of considerable interest is how the feedback yield depends on metallicity. Because star formation is expected to vary inversely with feedback yield, this has ramifications for predictions of SFRs in low-metallicity dwarfs in the local universe, as well as low-metallicity galaxies at high redshift. With the numerical implementation in TIGRESS-NCR, we can address this question; quantifying feedback yields over a wide range of metallicity and pressure is a key motivation for the present study.

Simulations that evolve the ISM with explicit treatments of physics and uniformly high resolution are valuable as laboratories where the interactions behind emergent properties of galaxies may be scrutinized in great detail. In addition, suites of such simulations offer a means to develop realistic subgrid treatments for deployment in large-scale, low-resolution galaxy formation models. While it is not possible to resolve the ISM directly in simulations of this kind, galactic-scale baryonic evolution depends entirely on the choices adopted for SFRs and galactic winds. An effort to develop a new generation of physically motivated, numerically calibrated subgrid models has been recently launched, under the umbrella of the SMAUG<fnref n="17" /> and Learning the Universe collaborations (see Kim et al. 2020a; Smith et al. 2024 as examples of subgrid wind modeling). In S. Hassan et al. (2024, submitted), a comparison was made between the SFR in galaxies in the Illustris-TNG50 simulation (Nelson et al. 2019; Pillepich et al. 2019), calculated with the Springel & Hernquist (2003) subgrid model, and what would have been predicted for the SFR using the PRFM theory and yield calibration from the TIGRESS-classic simulation suite (Ostriker & Kim 2022). This comparison shows intriguing differences: higher SFRs and shorter gas depletion times would be predicted from PRFM at high redshifts where gas is denser, and pressure is higher, compared to the native TNG SFR from the Springel & Hernquist (2003) subgrid model. However, for a fully quantitative and self-consistent prediction of the SFR (especially, at higher redshifts), it is critical to include a dependence of the yield on metallicity. Calibration of the feedback yield from TIGRESS-NCR simulations, which account for the dependence of both heating/cooling and radiative transfer on metallicity, is thus an important goal of this paper.

The remainder of this paper is organized as follows. We first summarize key methods in Section 2.1 and model parameters in Section 2.3. In Section 3, we then provide an overview for two chosen galactic conditions (solar neighborhood and conditions similar to the SFR-weighted mean in nearby star-forming galaxies from the PHANGS sample) with varying

<footnote n="16">¹⁶ Throughout this paper, we will use “metallicity” as a collective term for the total elemental abundance of heavy metals in gas and dust. When a distinction is needed, we will use “gas metallicity” and “dust abundance” to denote the metal abundance in the gas phase and dust, respectively.</footnote>

<footnote n="17">¹⁷ Simulating Multiscale Astrophysics to Understand Galaxies; https://www.simonsfoundation.org/flatiron/center-for-computational-astrophysics/galaxy-formation/smaug/.</footnote>
metallicities. We analyze the maps of gas and radiation properties and quantify the emergent SFRs as a function of metallicity. In Section 4, we analyze the full simulation suite in the context of the PRFM theory, and provide a new calibration to the feedback yield including metallicity dependence. We also introduce an effective equation of state for multiphase, star-forming gas in Section 4.4. We discuss and summarize our results in Sections 5 and 6.

## 2. Methods and Models

### 2.1. TIGRESS-NCR

We use the TIGRESS-NCR framework to run a suite of numerical simulations under widely varying galactic conditions, where the novel feature compared to our previous work is the range of gas metallicity and dust abundance. The TIGRESS-NCR framework models the star-forming ISM in a patch of a galactic disk including the effects of magnetic fields, gravity, galactic differential rotation, stellar feedback including UV radiation and SNe, cooling, heating, and chemistry.

We aim to treat most physical processes in the star-forming ISM as explicitly as possible (to the extent the resolution allows), avoiding ad hoc approximations. The numerical framework we have built has enabled many previous scientific studies, with technical details of the physical elements we have implemented described in several published papers. Interested readers should consult the TIGRESS-classic method paper (Kim & Ostriker 2017) for summaries of methods for MHD, shearing box, gravity, sink particles, and SN feedback (see also Kim et al. 2020b, for an update in gas accretion to sinks), and the TIGRESS-NCR method paper (Kim et al. 2023a) for an explanation of UV radiation transfer and review of selected photochemistry, cooling, and heating processes. Comprehensive descriptions of our formulations for photochemistry, cooling, and heating processes (separating neutral, photoionized, and collisionally ionized regimes), including detailed rate coefficients, are presented in Kim et al. (2023b), along with tests of our implementations.

Here, we briefly summarize the included chemical and thermodynamic processes that directly depend on gas metallicity and dust abundance.

1. Metal cooling. The metal cooling is directly proportional to gas metallicity. This includes fine-structure line cooling by C, C$^+$, O, rotational line cooling by CO, combined nebular line cooling in the warm ionized gas (Wolfire et al. 1995; Draine 2011; Vandenbroucke & Wood 2018), and metal cooling in the hot gas assuming collisional ionization equilibrium (Sutherland & Dopita 1993; Wiersma et al. 2009; Gnat & Ferland 2012).

2. Photoelectric heating on small grains and PAHs. Incident FUV photons cause electrons to be dislodged via the PE effect from surfaces of small grains and PAHs, sharing the excess energy with the surrounding gas (Watson 1972). The PE heating rate per hydrogen ($\Gamma_{\text{PE}}$) thus scales directly with the abundance of small grains, which we take as proportional to the total dust abundance ($Z'_d$); we do not explicitly consider possible variations of the PAH fraction (but see Draine et al. 2007; Aniano et al. 2020, for observational evidence of decreasing PAH fraction at low metallicities). The PE heating efficiency ($\epsilon_{\text{PE}}$) depends on grain charging, which depends on both the local radiation field ($J_{\text{FUV}}$) and the CR ionization rate ($\xi_{\text{cr}}$) via the free electron abundance $x_e$ (Bakes &

Tielens 1994; Wolfire et al. 1995; Weingartner & Draine 2001a).

3. Grain-assisted processes. We include grain-assisted recombination and H$_2$ formation that scale with dust abundance (Hollenbach & McKee 1979; Weingartner & Draine 2001a). The former also contributes to cooling. Heating and cooling related to H$_2$ formation and dissociation are included but play a minor role in overall energetics.

In our ART module, UV radiation in all three bands—PE (6 eV < $h\nu$ < 10.2 eV), Lyman–Werner (LW; 10.2 eV < $h\nu$ < 13.6 eV), and Lyman continuum (LyC; 13.6 eV < $h\nu$) —are attenuated by dust along each ray. In our adopted nomenclature, FUV means both PE and LW bands while EUV means the LyC band. We linearly scale the spectrum-averaged dust cross sections in Appendix B of Kim et al. (2023b) with dust abundance.

We also include heating by CR ionization that is independent of the metallicity. We adopt a simple prescription for the CR ionization rate $\xi_{\text{cr}} \propto \Sigma_{\text{SFR}}/\Sigma_{\text{gas}}$. The linear scaling with $\Sigma_{\text{SFR}}$ is motivated by the source of CRs in SNR shocks, while the inverse dependence on gas column density nominally represents collisional losses on large scales (as adopted by Wolfire et al. 2003). Strictly speaking, losses by transport out of the galactic disk likely exceed collisional losses at moderate ISM surface density (similar to the solar neighborhood), while following the above relation at higher surface density (most of our parameter space for the current simulation suite). We normalize based on solar neighborhood conditions, i.e., $\xi_{\text{cr},0} = 2 \times 10^{-16} \text{s}^{-1}$ is adopted for $\Sigma_{\text{SFR}} = 2.5 \times 10^{-3} \text{M}_\odot \text{kpc}^{-2} \text{yr}^{-1}$ and $\Sigma_{\text{gas}} = 10 \text{M}_\odot \text{pc}^{-2}$ (Indriolo et al. 2015). We apply an additional local attenuation recipe $\xi_{\text{cr}} \propto N_{\text{eff}}^{-1}$ if a local column density estimator $N_{\text{eff}} = 1.5 \times 10^{21} \text{cm}^{-2} (n_{\text{H}}/100 \text{cm}^{-3})^{0.3}$ exceeds $N_0 = 9.35 \times 10^{20} \text{cm}^{-2}$ (Neufeld & Wolfire 2017). We note that this column density estimator is different from what we adopted in Kim et al. (2023b) based on the ratio between the attenuated and unattenuated PE radiation fields from the radiation transfer solutions. The new form is adopted for simplicity and robustness since the previous prescription using the radiation transfer solutions does not converge with the ray truncation parameters.

The ART method we use in the TIGRESS-NCR framework is a direct but expensive method to follow the UV radiation fields. Even with highly optimized performance and parallel efficiency (Kim et al. 2017b), there are a few additional assumptions we adopt to reduce the overall cost: sparse ART calculation at every hydro time step of the warm and cold gas, and ray termination for FUV at $|z_{\text{p-p}}| = 300 \text{pc}$ above/below which we transition to a horizontally uniform FUV field from a plane-parallel radiation transfer solution<fnref n="18" /> (Kim et al. 2023a). In addition, two more parameters determine the ray termination conditions: (1) the maximum travel distance in the horizontal direction ($d_{\text{xy,max}}$) and (2) the ratio of the luminosity of the photon packet to the total luminosity of all sources in the domain ($\epsilon_{\text{PP}}$; FUV only). As we showed in Appendix of Kim et al. (2023a), the impact of these ray termination parameters on pressures, SFR surface density, and feedback yields are minimal, although the radiation fields at large distances from

<footnote n="18">¹⁸ There can be radiation source particles located above/below $z_{\text{p-p}}$. In such (very rare) cases, we neglect their FUV photons propagating upwards/downwards, but rays propagating downwards/upwards are followed consistently until the FUV photons are collected at the other side of $z_{\text{p-p}}$.</footnote>
<table><caption>**Table 1**
Input Physical Parameters</caption><thead><tr><th>Model Series<br>(1)</th><th>Z'<sub>d</sub><br>(2)</th><th>Σ<sub>gas,0</sub><br>(M<sub>☉</sub> pc<sup>-2</sup>)</th><th>Σ<sub>star</sub><br>(M<sub>☉</sub> pc<sup>-3</sup>)</th><th>z<sub>*</sub><br>(pc)</th><th>ρ<sub>dm</sub><br>(M<sub>☉</sub> pc<sup>-3</sup>)</th><th>Ω<br>(km s<sup>-1</sup> kpc<sup>-1</sup>)</th><th>R<sub>0</sub><br>(kpc)</th><th>L<sub>x,y</sub><br>(pc)</th><th>L<sub>z</sub><br>(pc)</th><th>Δx<br>(pc)</th></tr></thead><tbody><tr><td>S05</td><td>1, 0.1</td><td>5</td><td>1</td><td>500</td><td>0.002</td><td>15</td><td>8</td><td>2048</td><td>6144</td><td>8</td></tr><tr><td>R8</td><td>3, 1, 0.3, 0.1, 0.025</td><td>12</td><td>42</td><td>245</td><td>0.0064</td><td>28</td><td>8</td><td>1024</td><td>6144</td><td>8</td></tr><tr><td>S30</td><td>1, 0.1</td><td>30</td><td>42</td><td>245</td><td>0.0064</td><td>28</td><td>8</td><td>1024</td><td>6144</td><td>8</td></tr><tr><td>LGR4</td><td>3, 1, 0.3, 0.1, 0.025</td><td>50</td><td>50</td><td>500</td><td>0.005</td><td>60</td><td>4</td><td>512</td><td>3072</td><td>4</td></tr><tr><td>S100</td><td>1, 0.1</td><td>100</td><td>50</td><td>500</td><td>0.005</td><td>60</td><td>4</td><td>512</td><td>3072</td><td>4</td></tr><tr><td>S150</td><td>1, 0.1</td><td>150</td><td>50</td><td>500</td><td>0.005</td><td>100, 200</td><td>2</td><td>512</td><td>3072</td><td>4</td></tr></tbody></table>

**Note.** We assume $Z'_g = Z'_d$ except for $Z'_d = 0.025$ for which we adopt $Z'_g = 0.1$. In each model series, a suffix in the model name is used to denote the metallicity parameters, i.e., “-ZXX” means $Z'_g = Z'_d = XX$, or “-ZgXX, ZdYY” stands for $Z'_g = XX$ and $Z'_d = YY$. For R8 and LGR4, the additional suffix “b10” is used to denote models with weaker initial magnetic fields (initial plasma beta $\beta_0 = 10$ instead of $\beta_0 = 1$). For S150, an additional suffix for the galactic rotation parameters is added; “Om100q0” and “Om200” for $\Omega = 100 \text{ km s}^{-1} \text{ kpc}^{-1}$ with $q = 0.01$ and $\Omega = 200 \text{ km s}^{-1} \text{ kpc}^{-1}$, respectively. See Tables 2 and 3 for the expanded model names. The solar metallicity R8 and LGR4 models (R8-Z1 and LGR4-Z1) are identical to the R8 and LGR4 models presented in Kim et al. (2023a).

the midplane are more sensitive to these parameters. To save the computational cost, therefore, we first run the early stage of simulations (typically <1-2$t_{orb}$) using smaller $d_{xy,max}$ and larger $\epsilon_{PP}$ parameters. We then restart the simulation for a longer period (typically >2-4$t_{orb}$) with the values that give reasonable convergence for radiation fields in each model.

## 2.2. Simulation Initialization

For initial conditions, we adopt horizontally uniform, vertically stratified gas profiles for density, temperature, and magnetic field strength following double Gaussians representing warm and hot components (Kim & Ostriker 2017; Kim et al. 2023a). Following our previous practice, we also apply initial velocity perturbations with amplitude of $10-30 \text{ km s}^{-1}$ (higher values for higher surface density models) and create initial star clusters that provide initial UV radiation and SNe during the early evolution, prior to the formation of the first generation star clusters. The magnetic field is initialized along the local azimuthal direction ($\hat{y}$) with a vertically constant ratio of the initial magnetic pressure to the thermal pressure $\beta_0 \equiv 8\pi P_{th}/B_0^2$.

Galactic rotation is modeled in the local shearing box approximation (Stone & Gardiner 2010). The simulation box is centered on galactocentric radius $R_0$ and corotates with this point in the disk at an angular frequency $\Omega(R_0) = 2\pi/t_{orb}$. Galactic differential rotation is characterized by the shear parameter $q \equiv -d \ln \Omega/d \ln R|_{R_0}$. The resulting background flow along the local azimuthal direction ($\hat{y}$) is $v_0 = -q\Omega x\hat{y}$, sheared in the local radial coordinate (x); this flow is imposed in the initial conditions and maintained through shearing periodic boundary conditions at the x-faces of the box. The additional forces arising in this rotating frame are the Coriolis force $-2\Omega \times v$ and the tidal force $2q\Omega^2 x\hat{x}$. We assume a flat rotation curve ($q=1$) unless stated otherwise.

As we shall show, the final self-regulated state is insensitive to the initial conditions in most simulations. The initial gas profiles matter for the initial transient phase that typically lasts less than one orbit time ($t_{orb}$). Still, there remains a longer-term memory of initial magnetic fields (e.g., Kim & Ostriker 2015b), which converge to a self-consistent saturation state after a few orbit times unless the chosen initial field strength is widely different from the converged value. We defer to future investigations of the effects of initial magnetic fields and their growth due to galactic dynamo at varying rotation parameters and feedback.

## 2.3. Model Parameters

In the current simulation suite, the main gas parameters that matter the most are the initial gas surface density $\Sigma_{gas,0}$, gas metallicity $Z'_g$, and dust abundance $Z'_d$. Here, the prime means that values for the metal and dust ratios to gas mass are defined relative to solar neighborhood values; we adopt metal-to-gas $Z_{g,\odot} = 0.014$ (Asplund et al. 2009) and dust-to-gas $Z_{d,\odot} = 0.0081$ (Weingartner & Draine 2001b).

We assume $Z'_g = Z'_d$ for $Z'_g = 3, 1, 0.3$, and $0.1$ except for one model with $Z'_g = 0.1$ and $Z'_d = 0.025$. This choice is motivated by the observational evidence (e.g., Rémy-Ruyer et al. 2014; De Vis et al. 2019; Roman-Duval et al. 2022) showing the linear relation between $Z'_g$ and $Z'_d$ (or a constant dust-to-metal ratio) at near solar metallicities and a significant drop of the dust-to-metal ratio below $Z'_g = 0.1$. Although we trace the gas metallicity field locally using passive scalars with additional metal injection from SN ejecta, we do not use this cell-by-cell gas metallicity information for cooling and heating. We instead simply assume globally constant $Z'_g$ and $Z'_d$ for each simulation. Investigation of the effect of locally varying gas metallicity and dust abundance is deferred to future models with an explicit dust evolution (e.g., Hu et al. 2019, 2023b).

In addition to the self-gravity of gas, the vertical stratification of the gas disk is controlled by the vertical gravity of the old stellar disk and dark matter halo. We use a fixed potential for this “external” gravity adopting a functional form similar to Kuijken & Gilmore (1989) with three main parameters (see Kim & Ostriker 2017): stellar surface density $\Sigma_{star}$, stellar disk scale height $z_*$, and dark matter volume density $\rho_{dm}$.

We summarize the model parameters in Table 1. For the model name, we follow the naming convention used in Kim et al. (2023a) for galactic conditions with $\Sigma_{gas}$ and background gravity similar to solar neighborhood (R8), and similar to PHANGS galaxies (LGR4). For other models, we introduce a new naming convention simply representing initial gas surface density, e.g., S100 for $\Sigma_{gas,0} = 100 M_{\odot} \text{ pc}^{-2}$. For these models, $\Sigma_{gas}$ is similar to the TIGRESS-classic suite (Kim et al. 2020b; Ostriker & Kim 2022), but external gravity parameters are different.

Galactic rotation is parameterized by $\Omega$ and $q$. For most model series, we adopt a single value of $\Omega$, increasing with $\Sigma_{gas}$. We only vary the galactic rotation parameters for high gas surface density models where they make notable differences.

The parameters for different model series are chosen to cover a range of conditions in typical star-forming galactic disks. We then
explore the full metallicity variation for R8 and LGR4 models and only run two metallicities for other models. We use a single suffix with “Z” to denote the metallicity parameter models if $Z'_g = Z'_d$. The models with $Z'_g = 0.1$ and $Z'_d = 0.025$ use $Zg0.1Zd0.025$ as a suffix. Finally, we also explore the effect of initial magnetic field strength by running the R8 and LGR4 models with weaker initial fields (“b10” to denote $\beta_0=10$ while $\beta_0=1$ in other models) and the effect of galactic rotation in the S150 models (“Om100q0” for $\Omega=100$ km s$^{-1}$ kpc and $q=0.01$, close to a near rigid body rotation). For example, R8-Z1.0 is identical to R8-8pc in Kim et al. (2023a), while R8-Zg0.1Zd0.025-b10 denotes the same model with low metallicity $Z'_g = 0.1$ and dust abundance $Z'_d = 0.025$ with a weaker initial magnetic field.

The simulation domain is a vertically extended cuboid with the horizontal dimension of $L_x=L_y=512$ to 2048 pc and the vertical dimension of $L_z=3072$ to 6144 pc (see parameters in Table 1). We adopt a cubic resolution element with the side length of $\Delta x=4$ or 8 pc, depending on the models. The resolution convergence is demonstrated in Kim et al. (2023a) for R8 and LGR4 with $Z' = 1$.

# 3. Simulation Overview

## 3.1. Metallicity Dependence of the ISM Phases and Radiation Fields

Figures 1 and 2 qualitatively depict the effect of metallicity on the properties of gas and radiation fields for the R8 and LGR4 model series, respectively. From top to bottom, we present the maps from midplane ($z=0$) slices of hydrogen number density, gas phase, FUV (FUV = PE+LW; 6 eV $< h\nu < 13.6$ eV) radiation intensity (normalized relative to the mean intensity of the Draine field; Draine 1978, $J_{FUV,Draine} = 2.1 \times 10^{-4}$ erg s$^{-1}$ cm$^{-2}$ sr$^{-1}$), and radiation energy density of the LyC band. The gas metallicity and dust abundance decrease from left to right (from $(Z'_g, Z'_d) = (3, 3)$ to $(0.1, 0.025)$).

For the phase separation shown in the second row, we use a reduced definition from the original phase definition presented in Kim et al. (2023a). We use the temperature boundaries of $T=500$ K, 6000 K, 35000 K, and $5 \times 10^5$ K to divide the gas into the cold neutral medium (CNM), unstable neutral medium (UNM), warm medium (WNM+WIM), warm-hot ionized medium (WHIM), and hot ionized medium (HIM). For the warm medium, an additional hydrogen abundance cut is used to separate the warm neutral medium (WNM; $x_{H^0} > 0.5$) and warm ionized medium (WIM; $x_{H^0} < 0.5$).

The snapshots within each model series are chosen to have similar total luminosity emitted from sources, near a local star formation peak for each model. Yet, the FUV radiation field is much lower in the high-metallicity models (right two columns) with clear shadows cast by dense gas near sources. Only very young clusters ($t_{age} < 5$ Myr) are significant LyC sources. The LyC radiation field shows much sharper cutoffs due to the strong absorption by the neutral hydrogen. The fraction of LyC photons absorbed by dust also decreases as the dust abundance decreases allowing for more mass/volume to be ionized at a given ionizing luminosity (visually clearer in Figure 2). Although the gas density distribution is not very different, more pervasive UV radiation at lower metallicities allows more gas in the warm and hot phases than in the cold phase.

## 3.2. Star Formation Rates and Gas Depletion Times

Figure 3 shows the time evolution of SFR surface density for all members of the R8 (top) and LGR4 (bottom) model series, representative of galactic conditions with gas and stellar surface density comparable to the Milky Way and nearby disk galaxies. The SFR surface density is measured from the total star cluster mass with age less than $t_{bin}$, i.e.,

$$ \Sigma_{\text{SFR}} = \frac{\sum_i m_{*,i} (t_{\text{age}} < t_{\text{bin}})}{L_x L_y t_{\text{bin}}}. \qquad (1) $$

Here, we use $t_{bin}=40$ Myr. From these temporal histories, it is evident that SFRs are lower and also more bursty at lower metallicity.

The second and third columns compare $\Sigma_{SFR}$ and the gas depletion time

$$ t_{\text{dep}} \equiv \frac{\Sigma_{\text{gas}}}{\Sigma_{\text{SFR}}} \qquad (2) $$

as a function of metallicity (metallicity increases from left to right) during the late-stage evolution $t > 2t_{orb}$. It is immediately noticeable that SFR decreases as gas metallicity decreases, while, at $Z'_g = 0.1$, the further reduction of dust abundance does not make a difference in the SFR. Since long-term evolution with different SFR surface densities leads to different gas surface densities (see Table 2), the gas depletion time provides a better quantitative comparison, with $t_{dep}$ increasing from $\sim 2$ Gyr at $Z'_g = 3$ to $\sim 10$ Gyr at $Z'_g = 0.1$ in the R8 series, and a factor $\sim 3-4$ lower $t_{dep}$ in the LGR4 series.

# 4. PRFM with Varying Metallicities

To understand the metallicity dependence of SFRs seen in Section 3.2, we analyze the simulation results in the context of the PRFM theory of the star-forming ISM. We begin by demonstrating that in all cases the system is in vertical dynamical equilibrium, meaning that the total midplane pressure provides the support required to balance the total gravitational weight. We show also that the gas weight is insensitive to changes in metallicity (Section 4.2). This implies that the demand for pressure support remains the same, irrespective of metallicity. We then show that the thermal pressure (and magnetic stress to a lesser extent) increases at lower metallicity for a given SFR, meaning that low-metallicity systems are more efficient at recovering energy losses via stellar feedback (Section 4.3). We quantify the metallicity dependence of the feedback yield for each pressure/stress component. Finally, we provide an effective equation of state, the relation between the gas density and total pressure, averaged over the simulation domain (Section 4.4). Using the effective equation of state, we also compute effective vertical velocity dispersions; these are insensitive to metallicity, which is part of the reason why the weight is metallicity insensitive. The PRFM analysis explains the physics underpinning longer gas depletion times measured at low metallicity (Figure 3): a higher yield in the conversion of stellar feedback energy to ISM pressure at low $Z'_g$ and $Z'_d$ means that energy lost to dissipation in the ISM can be recouped at a lower SFR.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Visualization of the metallicity dependence of gas and radiation properties in the R8 model series (with $\Sigma_{gas,0} = 12 M_{\odot} \text{pc}^{-2}$ and $\Sigma_{star} = 42 M_{\odot} \text{pc}^{-2}$ comparable to the solar neighborhood); $Z'_g$ and $Z'_d$ increase from left to right. From top to bottom, we show midplane slices of number density, gas phase, FUV mean intensity normalized to the Draine field, and EUV radiation energy density. Young star clusters (color coded from age = 0 (magenta) to 40 Myr (cyan)) within $|z| < 50$ pc are overplotted as circles in the top row. We select snapshots from times near local star formation peaks with similar total luminosities from sources. From left to right, the snapshots are at $t = 453, 378, 420, 438, \text{ and } 428$ Myr with $\Sigma_{SFR} = 6.1, 3.4, 4.0, 3.0, \text{ and } 4.2 \times 10^{-3} M_{\odot} \text{ kpc}^{-2} \text{ yr}^{-1}$, $L_{LyC}/(L_{Ly}) = 6.3, 5.4, 7.5, 6.6, \text{ and } 6.5 L_{\odot} \text{ pc}^{-2}$, and $L_{FUV}/(L_{Ly}) = 19.2, 12.5, 13.7, 13.2, \text{ and } 15.5 L_{\odot} \text{ pc}^{-2}$.</figcaption></figure>

## 4.1. Definition of Measured Quantities

We first define and measure a variety of quantities that describe the vertical momentum conservation (i.e., vertical dynamical equilibrium in a steady state) and characterize the ISM disk. We construct horizontally averaged vertical profiles. We use the angle brackets $\langle q \rangle \equiv \int qdx dy/A$ with $A = L_x L_y$ to denote the horizontal average of any quantity $q$. To separately measure quantities in each phase, we use the Heaviside step function $\Theta(\text{ph})$ that returns 1 for the cell satisfying each phase definition or 0 otherwise (see Section 3.1 for the phase definition). For example, the area filling factor of each phase at a given height $z$ can be defined as

$$f_{A,ph}(z) = \langle \Theta(\text{ph}) \rangle. \qquad (3)$$

We then define the vertical profiles of “typical” density and pressure components, respectively, by

$$\rho_{\text{ph}}(z) = \frac{\langle \rho \Theta(\text{ph}) \rangle}{f_{A,ph}(z)} \qquad (4)$$

and

$$P_{\text{comp, ph}}(z) = \frac{\langle P_{\text{comp}} \Theta(\text{ph}) \rangle}{f_{A,ph}(z)}, \qquad (5)$$

where comp = th, turb, and mag to respectively denote the thermal component $P_{\text{th}}$, turbulent component (Reynolds stress) $P_{\text{turb}} \equiv \rho v_z^2$, and magnetic component (vertical Maxwell stress) $\Pi_{\text{mag}} \equiv B^2/(8\pi) - B_z^2/(4\pi)$ (e.g., Boulares & Cox 1990). In addition to the total Maxwell stress $\Pi_{\text{mag}}$, we also measure a decomposition of the stress into the mean $\Pi_B$ and turbulent $\Pi_{\delta B}$ using $\bar{B} \equiv \langle \mathbf{B} \rangle$ and $\delta \mathbf{B} \equiv \mathbf{B} - \bar{\mathbf{B}}$. The total pressure is then simply

$$P_{\text{tot,ph}} = P_{\text{th,ph}} + P_{\text{turb,ph}} + \Pi_{\text{mag,ph}}. \qquad (6)$$

Next, we calculate the weight profile $W(|z|) = (W_+ + W_-)/2$ using the integrals from a distance $|z|$ away from the midplane to
<figure><img src="image_5.png" /><figcaption>Figure 2. Same as Figure 1 but for LGR4 model series (with $\Sigma_{gas,0} = 50 M_{\odot} pc^{-2}$ and $\Sigma_{star} = 50 M_{\odot} pc^{-2}$, comparable to mean PHANGS galaxies). From left to right, the snapshots are at $t = 386, 233, 293, 333$, and 277 Myr with $\Sigma_{SFR} = 5.7, 5.5, 5.4, 5.8$, and $4.3 \times 10^{-2} M_{\odot} kpc^{-2} yr^{-1}$, $L_{Lyc}/(L_{Ly}) = 110, 95, 97, 79$, and $69 L_{\odot} pc^{-2}$, and $L_{FUV}/(L_{Ly}) = 272, 228, 242, 192$, and $191 L_{\odot} pc^{-2}$.</figcaption></figure>

the top and bottom of the box:

$$W_{\pm}(|z|) = \int_{\pm|z|}^{\pm L_z/2} \langle \rho \frac{d\Phi}{dz} \rangle dz, \quad (7)$$

where $\Phi$ is the total gravitational potential, which consists of the gravitational potential of gas, stars (young and old), and dark matter. We often decompose the weight into two terms; the weight from gas self-gravity ($W_{sg}$) and from the external gravity of stars and dark matter ($W_{ext}$). Where $W$ is used without an argument, it denotes the midplane ($z=0$) value.

Additionally, we measure the mass-weighted mean effective vertical velocity dispersion and gas scale height of the warm-cold two-phase gas ($2p=CNM+UNM+WNM$) defined by, respectively,

$$\sigma_{z,\text{eff}} \equiv \left( \frac{\int \langle P_{\text{tot}} \Theta(2p) \rangle dz}{\int \langle \rho \Theta(2p) \rangle dz} \right)^{1/2} \quad (8)$$

and

$$H = \left( \frac{\int \langle \rho \Theta(2p) \rangle z^2 dz}{\int \langle \rho \Theta(2p) \rangle dz} \right)^{1/2}. \quad (9)$$

Note that the effective velocity dispersion is a quadratic sum of thermal ($\sigma_{th}$), turbulent ($\sigma_{z,turb}$), and magnetic ($\sigma_{z,mag}$) components. Each of these may be defined by using Equation (8) but with $P_{th}$, $P_{turb}$, and $\Pi_{mag}$ instead of $P_{tot}$.

We summarize quantities measured from our simulations in Tables 2 and 3; values we report are averaged over a period during which the evolution is in a quasisteady state. Column (2) of Table 2 lists time ranges over which the time averages are taken, in units of Myr. We use a snapshot interval of ~1 Myr for the S05, R8, and S30 models and ~0.5 Myr for the LGR4, S100, and S150 models, providing typically ~500 snapshots for each model (~300 for S150 and 1000 for S05). Columns (3)–(5) list the mean gas surface density, SFR surface density, and gas depletion time discussed in Section 3.2, now including all models. Column (6) lists the mean star particle surface density. Columns (7)–(12) list the midplane hydrogen
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Summary of the emergent SFRs in the (a) R8 (top row) and (b) LGR4 (bottom row) model series, which have gas surface density comparable to local disk galaxies, showing lower SFRs at lower metallicities. Left: time evolution of Σ<sub>SFR</sub> from star cluster particles with age younger than 40 Myr. Early evolution for t < 2t<sub>orb</sub> (shaded area) is excluded in the analysis of the saturated state. Middle: box-and-whisker plots summarizing Σ<sub>SFR</sub> for the late stage evolution (t > 2t<sub>orb</sub>) as a function of Z'<sub>d</sub>. For visual clarity, the b10 models with initially weaker magnetic fields are shifted by 30% to the right. Right: box-and-whisker plots summarizing t<sub>dep</sub> ≡ Σ<sub>gas</sub>/Σ<sub>SFR</sub> for the late stage evolution (t > 2t<sub>orb</sub>) as a function of Z'<sub>d</sub>. Note that Σ<sub>gas</sub> at late stages are lower in higher metallicity models due to the prolonged high SFRs (see Table 2), making the difference in t<sub>dep</sub> larger than in Σ<sub>SFR</sub>. Initially, reduced magnetic fields (shown in the box-and-whisker plots with dashed lines) do not affect the late time evolution as all models reach similar saturated magnetic field strengths. At Z'<sub>d</sub> = 0.1, further reduction in the dust abundance (Z'<sub>d</sub> = 0.1 vs. Z'<sub>d</sub> = 0.025) does not make a difference in Σ<sub>SFR</sub>.</figcaption></figure>

number density, gas scale height, and effective velocity dispersions of the warm-cold two-phase gas. Note that the number density of hydrogen at the midplane is $n_{\text{mid,2p}} = \rho_{2p}(0)/(\mu_H m_H)$, where $\mu_H = 1.4$ is the mean molecular weight per hydrogen nucleus. Table 3 presents all midplane pressure and weight measurements that are used in the PRFM analysis in the following sections.

Figure 4 shows example midplane pressure/stress components of two selected metallicity models ($Z' = 0.1$ and 1) from the R8 and LGR4 series shown in Figures 1 and 2. Overall, the spatial variation of pressure/stress (in each component and total) within the two-phase gas is about an order of magnitude. Due to the diverging velocity field centered near the midplane, the hot superbubbles (bounded by the green contours) often have low turbulent stress. The magnetic stress also drops significantly within the hot gas—it can become negative if the vertical component dominates the magnetic field. The thermal pressure compensates somewhat to make the total pressure more uniform than each pressure/stress component even across the different gas phases. We note that these particular snapshots have $P_{\text{turb}}$ lower than is typical for the $Z' = 1$ models, where $P_{\text{turb}} \sim 2P_{\text{th}}$ on average (see Table 3). As we shall show in the subsequent sections, despite the large spatial and temporal

variation of pressure (e.g., Ostriker & Kim 2022; Kim et al. 2023a), a well-defined mean midplane pressure can be robustly determined, allowing us to investigate the quasisteady equilibrium state of the pressure/stress components.

A commonly adopted analytic weight estimator is

$$P_{\text{DE}} = \frac{1}{2} \pi G \Sigma_{\text{gas}}^2 + \Sigma_{\text{gas}} \sigma_{z,\text{eff}} (2G \rho_{\text{sd}})^{1/2}, \quad (10)$$

where the midplane stellar+dark matter volume density $\rho_{\text{sd}} = \Sigma_{\text{star}}/(2z_{*}) + \rho_{\text{dm}}$. In the last column in Table 3, we present results for $P_{\text{DE}}$, but take $\Sigma_{\text{gas}}^2 \rightarrow \Sigma_{\text{gas}}(\Sigma_{\text{gas}} + \Sigma_{\text{sp}})$ in the first term, where $\Sigma_{\text{sp}}$ is the surface density of star cluster particles formed in the simulation, since gas and young star clusters have similar vertical distributions. The contribution from $\Sigma_{\text{sp}}$ is only significant in the S100 and S150 models (see Table 2, column (6)). This weight estimator has been frequently adopted in observational work; an implicit assumption is that the gas disk's thickness is smaller than that of the stellar disk and dark matter halo.

For succinctness, in Tables 2 and 3 as well as the rest of the paper, we omit the subscript “2P” for the two-phase gas and the function argument “(0)” for the midplane value, i.e.,
<table><caption>**Table 2**

Time-averaged Properties in Simulations</caption><thead><tr><th>Model<br>(1)</th><th>(t<sub>1</sub>, t<sub>2</sub>)<br>(Myr, Myr)<br>(2)</th><th>Σ<sub>gas</sub><br>(M<sub>☉</sub> pc<sup>-2</sup>)<br>(3)</th><th>Σ<sub>SFR</sub><br>(M<sub>☉</sub> kpc<sup>-2</sup> yr<sup>-1</sup>)<br>(4)</th><th>t<sub>dep</sub><br>(Gyr)<br>(5)</th><th>Σ<sub>sp</sub><br>(M<sub>☉</sub> pc<sup>-2</sup>)<br>(6)</th><th>n<sub>H</sub><br>(cm<sup>-3</sup>)<br>(7)</th><th>H<br>(pc)<br>(8)</th><th>σ<sub>z,eff</sub><br>(km s<sup>-1</sup>)<br>(9)</th><th>σ<sub>z,turb</sub><br>(km s<sup>-1</sup>)<br>(10)</th><th>σ<sub>th</sub><br>(km s<sup>-1</sup>)<br>(11)</th><th>σ<sub>z,mag</sub><br>(km s<sup>-1</sup>)<br>(12)</th></tr></thead><tbody><tr><td>S05-Z1.0</td><td>(614, 1703)</td><td>4.07</td><td>1.13 · 10<sup>-4</sup></td><td>48.9</td><td>0.242</td><td>0.135</td><td>626</td><td>9.46</td><td>6.22</td><td>4.96</td><td>5.11</td></tr><tr><td>S05-Z0.1</td><td>(614, 1830)</td><td>4.76</td><td>4.74 · 10<sup>-5</sup></td><td>189</td><td>0.0616</td><td>0.167</td><td>494</td><td>8.17</td><td>3.84</td><td>6.25</td><td>3.58</td></tr><tr><td>R8-Z3.0</td><td>(438, 977)</td><td>8.74</td><td>3.74 · 10<sup>-3</sup></td><td>2.51</td><td>3.48</td><td>1.22</td><td>219</td><td>12.4</td><td>7.64</td><td>5.09</td><td>8.32</td></tr><tr><td>R8-Z1.0</td><td>(438, 977)</td><td>9.22</td><td>2.97 · 10<sup>-3</sup></td><td>3.79</td><td>2.55</td><td>1.16</td><td>235</td><td>12.8</td><td>8.38</td><td>5.92</td><td>7.61</td></tr><tr><td>R8-Z0.3</td><td>(438, 977)</td><td>9.96</td><td>1.92 · 10<sup>-3</sup></td><td>5.66</td><td>1.67</td><td>1.08</td><td>203</td><td>12.5</td><td>7.17</td><td>6.48</td><td>7.89</td></tr><tr><td>R8-Z0.1</td><td>(438, 977)</td><td>10.6</td><td>1.29 · 10<sup>-3</sup></td><td>14.8</td><td>0.970</td><td>1.23</td><td>188</td><td>11.9</td><td>6.36</td><td>6.60</td><td>7.60</td></tr><tr><td>R8-Zg0.1Zd0.025</td><td>(438, 977)</td><td>10.3</td><td>1.32 · 10<sup>-3</sup></td><td>12.2</td><td>1.20</td><td>1.28</td><td>200</td><td>12.1</td><td>6.89</td><td>6.55</td><td>7.47</td></tr><tr><td>R8-b10-Z1.0</td><td>(438, 977)</td><td>8.81</td><td>2.60 · 10<sup>-3</sup></td><td>3.77</td><td>2.97</td><td>1.14</td><td>219</td><td>12.6</td><td>7.50</td><td>5.98</td><td>8.18</td></tr><tr><td>R8-b10-Z0.3</td><td>(438, 977)</td><td>9.84</td><td>1.90 · 10<sup>-3</sup></td><td>6.90</td><td>1.88</td><td>1.26</td><td>203</td><td>12.2</td><td>7.17</td><td>6.39</td><td>7.52</td></tr><tr><td>R8-b10-Z0.1</td><td>(438, 977)</td><td>10.4</td><td>1.29 · 10<sup>-3</sup></td><td>15.1</td><td>1.18</td><td>1.19</td><td>200</td><td>12.2</td><td>6.48</td><td>6.66</td><td>7.91</td></tr><tr><td>R8-b10-<br>Zg0.1Zd0.025</td><td>(438, 977)</td><td>10.3</td><td>1.29 · 10<sup>-3</sup></td><td>15.6</td><td>1.32</td><td>1.09</td><td>231</td><td>13.4</td><td>7.08</td><td>6.56</td><td>9.26</td></tr><tr><td>S30-Z1.0</td><td>(438, 977)</td><td>16.2</td><td>0.0114</td><td>1.53</td><td>13.1</td><td>2.62</td><td>248</td><td>13.8</td><td>9.39</td><td>6.22</td><td>8.00</td></tr><tr><td>S30-Z0.1</td><td>(438, 977)</td><td>22.5</td><td>7.92 · 10<sup>-3</sup></td><td>3.47</td><td>6.67</td><td>3.03</td><td>212</td><td>13.6</td><td>9.07</td><td>6.93</td><td>7.46</td></tr><tr><td>LGR4-Z3.0</td><td>(204, 444)</td><td>35.7</td><td>0.0503</td><td>0.726</td><td>16.4</td><td>8.46</td><td>201</td><td>13.8</td><td>9.02</td><td>4.31</td><td>9.49</td></tr><tr><td>LGR4-Z1.0</td><td>(204, 488)</td><td>36.1</td><td>0.0343</td><td>1.08</td><td>14.6</td><td>7.38</td><td>171</td><td>13.6</td><td>8.19</td><td>5.57</td><td>9.26</td></tr><tr><td>LGR4-Z0.3</td><td>(204, 488)</td><td>42.0</td><td>0.0252</td><td>1.80</td><td>8.19</td><td>6.70</td><td>175</td><td>14.0</td><td>8.58</td><td>6.26</td><td>9.11</td></tr><tr><td>LGR4-Z0.1</td><td>(204, 488)</td><td>44.6</td><td>0.0186</td><td>3.16</td><td>5.43</td><td>7.28</td><td>184</td><td>13.7</td><td>7.90</td><td>6.51</td><td>9.06</td></tr><tr><td>LGR4-<br>Zg0.1Zd0.025</td><td>(204, 488)</td><td>42.4</td><td>0.0202</td><td>2.21</td><td>7.14</td><td>6.45</td><td>196</td><td>13.7</td><td>8.02</td><td>6.35</td><td>9.05</td></tr><tr><td>LGR4-b10-Z1.0</td><td>(204, 488)</td><td>33.2</td><td>0.0278</td><td>1.29</td><td>15.1</td><td>6.25</td><td>191</td><td>13.8</td><td>9.22</td><td>5.52</td><td>8.61</td></tr><tr><td>LGR4-b10-Z0.1</td><td>(204, 488)</td><td>41.4</td><td>0.0169</td><td>3.10</td><td>7.79</td><td>5.82</td><td>177</td><td>13.0</td><td>7.97</td><td>6.53</td><td>8.01</td></tr><tr><td>S100-Z1.0r</td><td>(204, 487)</td><td>54.4</td><td>0.127</td><td>0.495</td><td>42.9</td><td>14.2</td><td>276</td><td>18.4</td><td>15.7</td><td>5.77</td><td>7.83</td></tr><tr><td>S100-Z1.0</td><td>(204, 425)</td><td>42.0</td><td>0.0860</td><td>0.769</td><td>50.6</td><td>10.7</td><td>321</td><td>20.9</td><td>19.0</td><td>5.65</td><td>6.64</td></tr><tr><td>S100-Z0.1</td><td>(204, 488)</td><td>62.4</td><td>0.0885</td><td>0.799</td><td>34.7</td><td>12.4</td><td>215</td><td>16.3</td><td>12.4</td><td>6.95</td><td>8.01</td></tr><tr><td>S150-<br>Om200-Z1.0r</td><td>(153, 336)</td><td>105</td><td>0.247</td><td>0.453</td><td>49.4</td><td>36.2</td><td>113</td><td>17.1</td><td>6.50</td><td>4.87</td><td>15.1</td></tr><tr><td>S150-Om200-Z1.0</td><td>(153, 397)</td><td>88.7</td><td>0.133</td><td>0.743</td><td>66.0</td><td>21.1</td><td>157</td><td>20.0</td><td>6.57</td><td>4.86</td><td>18.3</td></tr><tr><td>S150-Om200-Z0.1</td><td>(153, 412)</td><td>115</td><td>0.108</td><td>1.10</td><td>37.5</td><td>30.9</td><td>111</td><td>17.0</td><td>5.17</td><td>6.32</td><td>14.9</td></tr><tr><td>S150-<br>Om100q0-Z1.0</td><td>(146, 307)</td><td>69.9</td><td>0.278</td><td>0.302</td><td>74.3</td><td>26.2</td><td>277</td><td>22.2</td><td>21.3</td><td>5.74</td><td>3.02</td></tr><tr><td>S150-<br>Om100q0-Z0.1</td><td>(146, 307)</td><td>68.7</td><td>0.261</td><td>0.319</td><td>74.5</td><td>21.4</td><td>214</td><td>17.7</td><td>15.7</td><td>7.20</td><td>3.74</td></tr></tbody></table>

Note. Column (1): model name, following convention explained in Table 1 and Section 2.3. Two models with a tag “r” (i.e., S100-Z1.0r and S150-Om200-Z1.0r) are the $Z' = 1$ simulations restarted from at the end of the early evolution for the $Z' = 0.1$ models. R8-Z3.0 and LGR4-Z3.0 are restarted from R8-Z1.0 and LGR4-Z1.0 at $t = 200$ and 150 Myr, respectively. Column (2): time interval over which the mean values are calculated. Columns (3)–(5): gas and SFR surface densities and their ratio, as presented in Section 3.2. Column (6): surface density of star cluster particles formed during the simulation. Columns (7)–(12): warm-cold two-phase (2p) gas properties as defined by Section 4.1.

$P_{tot,2p}(0) \to P_{tot}$, unless they are necessary to clarify the meaning. The weight integral without a function argument denotes the total weight integrated all the way to the midplane $z=0$ from the top/bottom of the box.

## 4.2. Vertical Dynamical Equilibrium

To investigate the behavior of SFR at different metallicities using the PRFM theory, we first check the establishment of vertical dynamical equilibrium. We again use the R8 and LGR4 model series for in depth investigations. Figure 5 shows the weight profiles along with the total pressure profiles for the 2p, WIM, and hot phases. Within one or two gas scale heights, the total pressure profiles of all phases are in good agreement with

the weight profiles and with each other. Only the midplane total pressure of the hot phase is smaller than those of the other phases. This is because the hot gas near the midplane represents the centers of expanding bubbles whose radial velocity increases outward. The excess pressure at high-$z$ in the 2p phase is dominated by the turbulent component arising from warm fountain flows (see also Kim & Ostriker 2018; Kim et al. 2020b).

Since the external vertical gravity is fixed at a given galactic condition, the total weights (or weight profiles) vary due to the change in total gas surface density, giving rise to a slightly higher weight at lower metallicity (see Table 2). Modulo this effect, the shape of total pressure and density profiles is more or less similar within each model series irrespective of metallicity, characterized by similar scale height and velocity dispersions.
<table><caption>**Table 3**
Measured Weight and Midplane Pressures in Simulations</caption><thead><tr><th>Model<br>(1)</th><th>W<br>(2)</th><th>P<sub>tot</sub><br>(3)</th><th>P<sub>turb</sub><br>(4)</th><th>P<sub>th</sub><br>(5)</th><th>Π<sub>mag</sub><br>(6)</th><th>Π<sub>dB</sub><br>(7)</th><th>Π<sub>B̄</sub><br>(8)</th><th>P<sub>tot,hot</sub><br>(9)</th><th>P<sub>DE</sub><br>(10)</th></tr></thead><tbody><tr><td>S05-Z1.0</td><td>0.125</td><td>0.141</td><td>0.0587</td><td>0.0400</td><td>0.0422</td><td>0.0120</td><td>0.0302</td><td>0.138</td><td>0.151</td></tr><tr><td>S05-Z0.1</td><td>0.146</td><td>0.130</td><td>0.0153</td><td>0.0883</td><td>0.0260</td><td>9.87 · 10<sup>-3</sup></td><td>0.0162</td><td>0.160</td><td>0.172</td></tr><tr><td>R8-Z3.0</td><td>1.61</td><td>1.69</td><td>0.618</td><td>0.246</td><td>0.829</td><td>0.371</td><td>0.458</td><td>1.37</td><td>1.75</td></tr><tr><td>R8-Z1.0</td><td>1.80</td><td>2.06</td><td>0.835</td><td>0.439</td><td>0.786</td><td>0.322</td><td>0.463</td><td>1.57</td><td>1.91</td></tr><tr><td>R8-Z0.3</td><td>1.84</td><td>2.01</td><td>0.571</td><td>0.575</td><td>0.862</td><td>0.327</td><td>0.535</td><td>1.41</td><td>2.04</td></tr><tr><td>R8-Z0.1</td><td>1.85</td><td>1.88</td><td>0.336</td><td>0.696</td><td>0.850</td><td>0.238</td><td>0.612</td><td>1.51</td><td>2.12</td></tr><tr><td>R8-Zg0.1Zd0.025</td><td>1.84</td><td>1.85</td><td>0.322</td><td>0.653</td><td>0.876</td><td>0.271</td><td>0.605</td><td>1.68</td><td>2.07</td></tr><tr><td>R8-b10-Z1.0</td><td>1.67</td><td>1.77</td><td>0.558</td><td>0.402</td><td>0.807</td><td>0.333</td><td>0.474</td><td>1.35</td><td>1.79</td></tr><tr><td>R8-b10-Z0.3</td><td>1.80</td><td>1.85</td><td>0.465</td><td>0.589</td><td>0.791</td><td>0.303</td><td>0.487</td><td>1.42</td><td>1.98</td></tr><tr><td>R8-b10-Z0.1</td><td>1.87</td><td>1.87</td><td>0.308</td><td>0.682</td><td>0.882</td><td>0.247</td><td>0.635</td><td>1.46</td><td>2.11</td></tr><tr><td>R8-b10-Zg0.1Zd0.025</td><td>1.96</td><td>1.94</td><td>0.324</td><td>0.597</td><td>1.02</td><td>0.268</td><td>0.752</td><td>1.57</td><td>2.24</td></tr><tr><td>S30-Z1.0</td><td>4.30</td><td>4.69</td><td>1.86</td><td>1.07</td><td>1.76</td><td>0.948</td><td>0.809</td><td>3.64</td><td>3.97</td></tr><tr><td>S30-Z0.1</td><td>5.67</td><td>6.21</td><td>2.03</td><td>2.04</td><td>2.14</td><td>0.959</td><td>1.18</td><td>4.46</td><td>5.91</td></tr><tr><td>LGR4-Z3.0</td><td>10.9</td><td>12.7</td><td>6.41</td><td>1.57</td><td>4.69</td><td>2.20</td><td>2.50</td><td>10.2</td><td>9.82</td></tr><tr><td>LGR4-Z1.0</td><td>10.4</td><td>11.7</td><td>5.00</td><td>2.43</td><td>4.31</td><td>1.86</td><td>2.45</td><td>10.2</td><td>9.56</td></tr><tr><td>LGR4-Z0.3</td><td>11.8</td><td>12.7</td><td>3.88</td><td>3.18</td><td>5.65</td><td>1.83</td><td>3.82</td><td>11.9</td><td>12.1</td></tr><tr><td>LGR4-Z0.1</td><td>12.6</td><td>13.2</td><td>3.37</td><td>3.97</td><td>5.91</td><td>2.07</td><td>3.83</td><td>11.4</td><td>13.1</td></tr><tr><td>LGR4-Zg0.1Zd0.025</td><td>11.8</td><td>13.3</td><td>5.00</td><td>3.24</td><td>5.07</td><td>1.89</td><td>3.18</td><td>10.1</td><td>12.1</td></tr><tr><td>LGR4-b10-Z1.0</td><td>8.99</td><td>10.1</td><td>4.23</td><td>2.00</td><td>3.82</td><td>1.34</td><td>2.48</td><td>9.04</td><td>8.55</td></tr><tr><td>LGR4-b10-Z0.1</td><td>10.9</td><td>11.5</td><td>3.44</td><td>3.41</td><td>4.62</td><td>1.66</td><td>2.96</td><td>10.2</td><td>11.5</td></tr><tr><td>S100-Z1.0r</td><td>27.7</td><td>32.8</td><td>19.6</td><td>5.70</td><td>7.43</td><td>3.20</td><td>4.23</td><td>24.1</td><td>21.0</td></tr><tr><td>S100-Z1.0</td><td>22.2</td><td>35.4</td><td>27.8</td><td>3.93</td><td>3.66</td><td>1.74</td><td>1.93</td><td>18.4</td><td>15.5</td></tr><tr><td>S100-Z0.1</td><td>28.2</td><td>31.6</td><td>12.7</td><td>8.53</td><td>10.4</td><td>5.24</td><td>5.13</td><td>24.9</td><td>24.1</td></tr><tr><td>S150-Om200-Z1.0r</td><td>70.9</td><td>72.2</td><td>15.9</td><td>8.81</td><td>47.5</td><td>14.3</td><td>33.2</td><td>69.1</td><td>56.3</td></tr><tr><td>S150-Om200-Z1.0</td><td>58.7</td><td>57.4</td><td>8.78</td><td>3.82</td><td>44.8</td><td>6.21</td><td>38.6</td><td>65.1</td><td>45.1</td></tr><tr><td>S150-Om200-Z0.1</td><td>72.7</td><td>72.0</td><td>8.66</td><td>12.4</td><td>51.0</td><td>9.55</td><td>41.4</td><td>70.1</td><td>64.7</td></tr><tr><td>S150-Om100q0-Z1.0</td><td>53.4</td><td>112</td><td>98.7</td><td>10.6</td><td>3.04</td><td>2.48</td><td>0.566</td><td>34.6</td><td>33.0</td></tr><tr><td>S150-Om100q0-Z0.1</td><td>42.0</td><td>53.1</td><td>32.0</td><td>16.2</td><td>4.86</td><td>3.79</td><td>1.08</td><td>39.5</td><td>29.1</td></tr></tbody></table>

Note. Columns (2)–(8): weight and pressure/stress for the warm-cold two-phase (2p) gas as defined in Section 4.1. Column (9): total pressure of the hotphase. Column (10): weight estimator (Equation (10)). All values are in units of $10^4 k_B \text{ cm}^{-3}$ K, where $k_B$ is the Boltzmann constant.

To give a more quantitative view, Figure 6(a) plots the measured total midplane pressure of the 2p phase as a function of the total weight for all models. This pressure and the weight are in excellent agreement, validating vertical dynamical equilibrium.<fnref n="19" /> For the same galactic condition, the models with different metallicities (denoted by different colors) are clustered at similar weights with a tendency to have slightly higher weights at lower metallicities. Again, this is because, at lower metallicity, the gas surface density decreases more slowly due to lower SFR surface density and longer gas depletion time.

Figure 6(b) shows that total pressure equilibrium is satisfied between the 2p and hot phases. Figure 6(c) shows that the analytic weight estimator in Equation (10) using the measured effective velocity dispersion (from the mass-weighed average along the vertical direction) provides a good estimate of the true weight and hence the total pressure. These conclusions

hold irrespective of metallicity. While above each panel we report the results of a bivariate fit of pressure to W and $Z'_g$, the purpose of this is simply to demonstrate (1) how accurately vertical equilibrium is satisfied, and (2) how negligible the measured dependence on metallicity is.

Taken together, and using the definition of weight in Equation (7), the results of Figure 6 show that, for all gas phases,

$$P_{\text{tot}} \approx \frac{1}{2} \sum_{\text{gas}} \langle g_z \rangle \quad (11)$$

where

$$\langle g_z \rangle \approx \pi G \sum_{\text{gas}} + 2\sigma_{z,\text{eff}} (2G\rho_{\text{sd}})^{1/2} \quad (12)$$

is an estimator of the mean vertical gravity within the gas layer. We caution, though, that the weight estimator in Equation (10) and gravity estimator in Equation (12) are derived assuming that the gas disk is thinner than the stellar disk and dark matter halo. This assumption, and hence application of Equation (10), is acceptable for the current simulation suite and most nearby normal star-forming galaxies (e.g., Sun et al. 2020; Barrera-Ballesteros et al. 2021). However, more extreme, gas-rich systems (e.g., Girard et al. 2021) may have more vertically

<footnote n="19">¹⁹ Only two models (S150-Om100q0-Z1.0 and S100-Z1.0) show significant deviation from equilibrium, with $P_{tot} > W$ by a factor of ~2.1 and 1.6, respectively. A deeper inspection of these models reveals that they are still in vertical dynamical equilibrium when we account for all phases, and the hot phase occupies the majority of the volume for all z. In this case, the warm/ cold 2p phase is no longer volume-filling and compressed into clumpy blobs, making the 2p phase not a representative phase for vertical dynamical equilibrium.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 4. Example pressure/stress component maps for two selected metallicity models ($Z' = 0.1$ and 1) for the R8 (left) and LGR4 (right) series using the same snapshot time as shown in Figures 1 and 2, respectively. The green contours denote the boundary between WIM and WHIM, equivalent to a temperature boundary at $T = 3.5 \times 10^4$ K.</figcaption></figure>

extended gas disks where Equation (10) overestimates the weight. For this reason, in the future, we recommend adopting the generalized formulae for weight and vertical gravity presented in S. Hassan et al. (2024, submitted), which may be applied to gas disks that are either thinner or thicker than the corresponding stellar disks.

Figure 7 shows the vertical profiles of pressure components in the 2p phase. Unlike the total pressure and weight profiles, these show noticeable metallicity dependence. In particular, the thermal pressure increases as the gas metallicity decreases, while there is a corresponding reduction in turbulent Reynolds stress. The Maxwell stresses have somewhat less consistent behavior; for R8, the mean Maxwell stress is higher, and the turbulent Maxwell stress is lower at low metallicity, while for LGR4 these trends are less clear.

We conclude that, consistent with expectations, the disk's vertical structure satisfies the quasisteady-state vertical momentum conservation law such that the total pressure support matches the weight, which is insensitive to metallicity. However, the detailed ISM physics governed by MHD and thermodynamics with radiation and mechanical feedback is responsible

for setting the individual pressure components. To fulfill the required total pressure support, the pressure components are constrained to add up to the same value regardless of metallicity, but the relative importance of different components varies in interesting ways. From now on, we shall quantify the variation of the relative importance of each pressure component as a function of weight and metallicity.

### 4.3. Feedback Yields

In the PRFM theory, we postulate that the ISM weight is balanced by pressure support that responds to star formation feedback, which can vary as needed to fulfill its role. The connection between pressure and SFR surface density can be characterized in terms of the feedback yield parameter (see Ostriker & Kim 2022, and references therein; and Kim & Ostriker 2015b for an in-depth focus on magnetic effects). The total feedback yield, and feedback yields of individual stress components, is defined by

$$ \Upsilon_{\text{tot}} \equiv \frac{P_{\text{tot}}}{\Sigma_{\text{SFR}}}; \quad \Upsilon_{\text{comp}} \equiv \frac{P_{\text{comp}}}{\Sigma_{\text{SFR}}} \qquad (13) $$
<figure><img src="image_8.png" /><figcaption>Figure 5. Demonstration of vertical dynamical equilibrium using the vertical profiles of the R8 (top row) and LGR4 (bottom row) models. Metallicity increases from left to right. We plot the vertical profiles of gas weight (W, thick solid), total pressure of the 2P (P<sub>tot,2P</sub>, thin solid), WIM (P<sub>tot,WIM</sub>, dotted), and hot (P<sub>tot,hot</sub>, dashed) phases. The lines and shaded areas show the median and 16th to 84th percentiles for the time range in Table 2. Note that, for visual clarity, we only show z = ± L<sub>z</sub>/4 rather than the full extent of z = ± L<sub>z</sub>/2.</figcaption></figure>

using the pressure/stress at the midplane reported in Table 3. The feedback yield has units of velocity.

Figure 8 shows each feedback yield measured from all models as a function of the weight and metallicity. We fit the resulting feedback yield using a bi-variate power-law model for the weight and gas metallicity. The fitting is carried out by the orthogonal distance regression method implemented in the scipy package. Since the goal here is to provide the calibration of the feedback yield employing the equilibrium assumption made in the PRFM theory,<fnref n="20" /> we simply use the mean values measured over long-term evolution covering four to six star formation/feedback and outflow/inflow cycles and weigh them equally. We do not attempt to reduce the bias arising from the uneven sampling of the parameter space given that the underlying true parameter distribution is not well understood. We remark that, if the true physical relationship follows this theoretical power law, then this regression should converge to the true parameters regardless of parameter priors. The bivariate fitting results are presented at the top of the following figures.

The fitting results are

$$ \Upsilon_{\text{th}} = 390 \text{ km s}^{-1} W_4^{-0.46} Z_g^{-0.53}, \quad (14) $$

$$ \Upsilon_{\text{turb}} = 561 \text{ km s}^{-1} W_4^{-0.21} Z_g^{-0.04}, \quad (15) $$

$$ \Upsilon_{\text{mag}} = 578 \text{ km s}^{-1} W_4^{-0.40} Z_g^{-0.44}, \quad (16) $$

where $W_4 \equiv W/(10^4 k_B \text{ cm}^{-3} \text{ K})$. In the high-pressure regime, the different galactic rotation speeds lead to diverging magnetic

and turbulent feedback yields (upper versus lower triangles). There is a rough correspondence between low(high) $\Upsilon_{\text{mag}}$ and high(low) $\Upsilon_{\text{turb}}$, implying an exchange between turbulent and magnetic energy. It is thus reasonable to combine turbulent and magnetic stresses and measure the associated feedback yield as in Figure 8(d), which gives

$$ \Upsilon_{\text{turb+mag}} = 1.17 \times 10^3 \text{ km s}^{-1} W_4^{-0.22} Z_g^{-0.18}. \quad (17) $$

The decreasing trend of $\Upsilon_{\text{th}}$ with increasing $W = P_{\text{tot}}$ and $Z'_g$ is the cleanest relation seen in Figure 8. As we will show below, this behavior is mainly driven by the increased attenuation of FUV radiation in environments with higher density (corresponding to higher weight and pressure) and higher dust abundance. A decreasing trend of $\Upsilon_{\text{th}}$ with increasing pressure and weight at solar metallicity was previously reported in Ostriker & Kim (2022) based on the TIGRESS-classic suite at $Z' = 1$ (shown as a dashed line in Figure 8(a)). In the TIGRESS-classic framework, however, the global radiation attenuation effect is imposed by a theoretical model for the radiation field based on a simple plane-parallel approximation (see Equation (11) in Ostriker & Kim 2022). With explicit UV radiation transfer using ART as implemented in TIGRESS-NCR, in Kim et al. (2023a), we confirmed for R8 and LGR4 (at solar metallicity) consistency in the scaling of $\Upsilon_{\text{th}}$ with $W$ with that of the TIGRESS-classic suite, with a normalization about a factor of 1.5 higher. Figure 8(a) reaffirms the evidence that $\Upsilon_{\text{th}}$ decreases with $W$, but now over a wider range of parameter space and at different metallicities. This is an indirect confirmation of the effective attenuation model used in the TIGRESS-classic framework (Kim et al. 2020b; Ostriker & Kim 2022); a direct comparison of the resulting FUV radiation field with different approximate solutions will be presented in N. Linzer et al. (2024, submitted).

<footnote n="20">²⁰ It is of great interest to develop a theoretical model including the effects of mutually correlated time evolution of pressure and SFR surface density. However, such a model requires a deeper understanding of the timescales of energy gain and loss as well as dynamical evolution under changing gravity and pressure, which we defer to future work.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Demonstration of vertical dynamical equilibrium for all models. We plot the measured total midplane pressures of the 2p and hot phases (a) $P_{tot}$ and (b) $P_{tot,hot}$, respectively, as well as (c) the weight estimator $P_{DE}$ presented in Equation (10) as a function of the total weight $W$. The bottom panel gives relative differences between the ordinates and abscissas. Each point represents the mean value of each model over the time range in Table 2. The color of points shows the gas metallicity with circles for $Z'_g = Z'_d$ and squares for $Z'_d = 0.025$. The star symbols are for the b10 models (they are hardly visible as they overlap with circles so well). The lower and upper triangles are for the s150-0m100q0 and s150-0m200 models. The colored lines are the fitting results of a bi-variate power-law model. The dashed black line is for the one-to-one relation.</figcaption></figure>

To understand the metallicity dependence of the thermal feedback yield, multiple factors must be considered simultaneously. The simplest estimation of the thermal pressure is from the assumption of the balance between metal cooling ($n_H \Lambda_{metal}$, especially C II and O I) and grain PE heating ($n_H \Gamma_{PE}$) in the CNM, and pressure equilibrium between the CNM and WNM. The assumption of thermal equilibrium gives rise to $n_{eq} = \Gamma_{PE}/\Lambda_{metal}$, and hence, one can write

$$P_{\text{th,eq}}/k_B = n_{\text{eq}}T = \frac{\Gamma_{\text{PE}}T}{\Lambda_{\text{metal}}(T, x_c)} \propto \frac{\epsilon_{\text{PE}}J_{\text{FUV}}Z'_d}{Z'_g} \quad (18)$$

where in the final proportionality we focus just on the factors that depend strongly on metallicity. Dividing Equation (18) by $\Sigma_{SFR}$, we obtain

$$\gamma_{\text{th}} \propto \frac{Z'_{d}}{Z'_{g}} \epsilon_{\text{PE}} f_{\tau}. \quad (19)$$

Here, $f_{\tau} \equiv 4\pi J_{FUV}/\Sigma_{FUV}$ is an attenuation factor relating the FUV radiation field to the injection rate of FUV per unit area in the disk, $\Sigma_{FUV}$, which is linearly proportional to $\Sigma_{SFR}$ averaged over an appropriate time bin (e.g., $t_{bin} \sim 10$ Myr for FUV).

Although the metallicity effects cancel out as long as $Z'_g$ varies linearly with $Z'_d$ to the lowest order, the higher order effects associated with radiation transfer ($f_{\tau}$) and the thermal and ionization state of gas ($\epsilon_{PE}$) still lead to a metallicity dependence of the thermal feedback yield in the CNM. On the one hand, as clearly seen in Figures 1 and 2, the FUV radiation field (third row; $\chi_{FUV} \equiv J_{FUV}/J_{FUV, Draine}$) is much higher at low dust abundance (with similar $\Sigma_{SFR}$) because attenuation is vastly reduced, making $f_{\tau}$ significantly higher. On the other hand, $\epsilon_{PE}$ depends (inversely) on the grain charging parameter $\propto \chi_{FUV} T^{1/2}/n_e$ (Bakes & Tielens 1994; Weingartner & Draine 2001a), with $\epsilon_{PE}$ larger in colder, higher density gas,

and sensitive to the ionization fraction. In the neutral ISM, the main source of free electrons is from H ionization by low-energy CRs (with electrons from $C^+$ becoming important at $Z'_g > 1$), so the grain charging parameter is relatively insensitive to metallicity for a given FUV field and CR ionization rate (see Figure 6 of Kim et al. 2023b). In R8 and LGR4, we find that the reduction from $Z'_g = 1$ to 0.1 gives rise to a factor of 2.6 and 3 increase in $f_{\tau}$ and 0.9 and 0.6 decrease in $\epsilon_{PE}$, respectively. This explains the larger thermal feedback yield at lower metallicity.

Finally, there is no difference seen in the model with reduced dust abundance $Z'_d = 0.025$ at $Z'_g = 0.1$ (square symbols). We find that this is because, at this low metallicity, the PE heating is comparable to or weaker than the CR heating, which has no dependence on dust abundance. At the same time, metallicity-dependent cooling is no longer dominant. Rather, the H I Lyα cooling in the WNM becomes dominant as the CNM fraction decreases with decreasing $Z'_g$, removing the metallicity dependence below $Z'_g < 0.1$.

The turbulent feedback yield shown in Figure 8(b) varies less sensitively with weight and metallicity. This is qualitatively consistent with the expectation that the total momentum injected by each SN is not a sensitive function of gas density and metallicity (e.g., Thornton et al. 1998; Kim & Ostriker 2015a; Martizzi et al. 2015; Kim et al. 2017a, 2023b; Karpov et al. 2020; Steinwandel et al. 2020; Oku et al. 2022). In a high-pressure environment, the behavior becomes somewhat irregular and more sensitive to other environmental parameters (e.g., galactic rotation speed). The apparent behavior depending on galactic rotation (upper and lower triangles at high pressure/weight) is related to magnetic fields.

In the models with high rotation/shear (S150-0m200; upper triangles), the mean magnetic field gets much stronger than in the low shear models (S150-0m100q0; lower triangles), dominating the overall support against the total weight ($\Upsilon_{mag} \sim 10^3 \text{ km s}^{-1}$)
<figure><img src="image_10.png" /><figcaption>Figure 7. Vertical profiles of pressure components in the 2p phase as defined by Equation (5). From left to right, we show turbulent, thermal, mean magnetic, and turbulent magnetic components, while different colors denote different metallicity in R8 (top) and LGR4 (bottom). The lines and shaded areas show the median and 16th to 84th percentiles for the time range in Table 2.</figcaption></figure>

with $\Upsilon_{\text{turb}}$ and $\Upsilon_{\text{th}} \sim 10^2 \text{ km s}^{-1}$). The gas depletion time is longer in S150-Om200, with a significantly lower turbulent velocity dispersion $\sigma_{z,\text{turb}} \sim 5-7 \text{ km s}^{-1}$ as opposed to 15–20 km s$^{-1}$ in S150-Om100q0. In the S100 and S150-Om100q0 models where SN feedback is very strong and blows away a lot of mass together with magnetic fields, the buildup of the mean magnetic field is hindered significantly. These models behave like hydrodynamics models without much contribution from the magnetic feedback yield ($\Upsilon_{\text{mag}} \lesssim 100 \text{ km s}^{-1}$), which is compensated by the large turbulent feedback yield.

The total feedback yield is the sum of all components. In Figure 9, we show the total feedback yield measurements with (a) a fitting result to the total feedback yield directly and (b) a two-component model using thermal and nonthermal feedback yields (i.e., the sum of Equations (14) and (17)). The direct fitting result is

$$ \Upsilon_{\text{tot}} = 1.65 \times 10^3 \text{ km s}^{-1} W_4^{-0.29} Z_g^{-0.27}. \quad (20) $$

There is only a small improvement in the mean of the residual from (b) compared to (a). This implies that these two are virtually equivalent within the parameter space explored in this paper. But, the two-component model would behave better in extreme conditions by respectively capturing the dominance of thermal and nonthermal pressure supports in low- and high-pressure regimes. Compared to the reference line showing the result from the TIGRESS-classic suite (Ostriker & Kim 2022), TIGRESS-NCR gives quite similar values near $W/k_B = 10^5 \text{ cm}^{-3} \text{ K}$ for $Z'_g = 1$ but has a steeper slope, with feedback yields larger by $\sim 0.2$ dex at low pressure. Additionally, $\Upsilon_{\text{tot}}$ is larger by $\sim 0.3$ dex at $Z'_g = 0.1$ compared to $Z'_g = 1$.

We attribute the enhancement of total feedback yield to the effect of ionizing UV, which was missing in TIGRESS-classic.

Further investigation of the role of individual feedback channels will be addressed in a subsequent paper.

Finally, Figure 10 shows the relationship between $\Sigma_{\text{SFR}}$, $W$, and $Z'_g$ for our full simulation suite, together with a bivariate fit.<fnref n="21" /> The superlinear dependence on $W$ and positive dependence on $Z'_g$ express the fact that all forms of feedback become less efficient in regions where the pressure (and density) and the metallicity are higher. Since, however, the dependencies on $W$ and $Z'_g$ in Equation (20) are comparable, while the variation in $W$ is larger than the variation in $Z'_g$ if one considers the range of conditions under which most of the stars in the universe have formed, the primary environmental factor affecting star formation is $W$.

### 4.4. Effective Equation of State

Figure 11 presents the mean relation between total midplane pressure and gas density of all gas (left) and just the 2p gas (right). The former relation represents an effective equation of state in a volume-averaged sense, encapsulating the system-level behavior of the multiphase star-forming ISM gas that develops as a consequence of the thermodynamic and MHD responses to intermittent energy injection by feedback. We find that the effective equation of state is not sensitive to changes in metallicity. The best-fit bi-variate power-law relation gives

$$ P_{\text{tot}}/k_B = 2.85 \times 10^4 \text{ cm}^{-3} \text{ K} \left( \frac{n_{\text{mid}}}{\text{cm}^{-3}} \right)^{1.29} Z_g^{-0.10}. \quad (21) $$

From the TIGRESS-classic simulations, we previously found a slightly steeper (1.43) slope when considering just the 2p gas; this comparison is shown in the right panel.

<footnote n="21">²¹ The fit is slightly different from what would be obtained from $\Sigma_{\text{SFR}} = W/\Upsilon_{\text{tot}}$ using Equation (20) because the fitting results are not precise at the level of the significant figure reported here.</footnote>
<figure><img src="image_11.png" /><figcaption>(a) Thermal Feedback Yield</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Turbulent Feedback Yield</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) Magnetic Feedback Yield</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) Non-Thermal Feedback Yield</figcaption></figure>

**Figure 8.** Metallicity and weight dependence of each feedback yield component. Symbols denote the mean value of feedback yield and weight over the time range in Table 2. The color of points shows the gas metallicity with circles for $Z'_g = Z'_d$ and squares for $Z'_d = 0.025$. The star symbols are for the b10 models. The lower and upper triangles are for the S150-Om100q0 and S150-Om200 models. The fitting results are presented as the equation in the panel title and as colored lines. The relative errors of the fitting results (i.e., model/simulation - 1) are presented at the bottom of each panel. In panels (a) and (b), the fits to the TIGRESS-classic suite (with $Z' = 1$) are shown in Ostriker & Kim (2022) as the black dashed line.
<figure><img src="image_15.png" /><figcaption>**Figure 9.** Two fitting results for the total feedback yield. Left: total feedback yield model obtained by direct fitting for measured $\Upsilon_{tot}$, i.e., Equation (20). Right: total feedback yield model obtained by the sum of two bivariate power fitting results in Figures 8(a) and (d), i.e., Equation (14) plus Equation (17). The fit to the TIGRESS-classic suite presented in Ostriker & Kim (2022) is shown as the black dashed line. Symbols and lines have the same meaning as in Figure 8.</figcaption></figure>

It should be borne in mind that the above effective equation of state represents a numerical average over large space and timescales, which physically corresponds to multiple cycles of star formation and feedback. Within a given local cycle, a smaller-scale region in a disk will experience first an increase in density and pressure (from compression driven by turbulence and gravity), and subsequently, an increase in specific energy (due to feedback) that leads to a drop in density. The effective equation of state represents an average over the complex dynamic and thermodynamic cycles that take place within the ISM under different galactic environmental conditions. By relating the quasiequilibrium midplane density and pressure values, the effective equation of state differs from an adiabatic pressure–density relation (applicable over time for a given fluid element) or polytropic pressure–density relation (applicable throughout space for a given equilibrium structure). The effective equation of state is more analogous to a polytropic relation, in that it can be used to compute the equilibrium thickness of the disk given its surface density, similar to the mass–radius relation for stellar polytropes.

The right panel of Figure 11 can be translated into an effective midplane vertical velocity dispersion of the 2p gas,

$$ \sigma_{z,\text{eff, mid}} = \left( \frac{P_{\text{tot, 2p}}}{\rho_{2p}} \right)^{1/2} = 8.9 \text{ km s}^{-1} \left( \frac{P_{\text{tot, 2p}}/k_B}{10^4 \text{ cm}^{-3} \text{ K}} \right)^{0.08} Z_g'^{-0.005}. \quad (22) $$

The midplane effective velocity dispersion can be compared with the mass-weighted mean considering the warm and cold

gas over the whole vertical domain, shown in the left panel of Figure 12, which has a power-law fit

$$ \sigma_{z,\text{eff}} = 11.7 \text{ km s}^{-1} W_4^{0.12} Z_g'^{0.03}. \quad (23) $$

The mass-weighted mean is ~30% higher than the midplane value. By comparison, the effective mass-weighted velocity dispersion reported from TIGRESS-classic simulations in Ostriker & Kim (2022) is comparable to Equation (23) for solar neighborhood conditions, but increases faster in high-pressure environments ($\propto W^{0.22}$). The lower effective velocity dispersion for TIGRESS-NCR in the high-weight regime is mainly because SN-driven turbulence ends up being reduced when ionizing radiation (early feedback) is included.

The right panel of Figure 12 shows the turbulent component of the mass-weighted mean velocity dispersion. Similar to the total, the turbulent component also scales weakly with weight and metallicity. In the high-pressure regime, there is a wide range of turbulent velocity dispersions, similar to the divergence among models of the turbulent feedback yield due to greatly enhanced magnetic support in rapidly rotating models.

### 4.5. Application to Subgrid Modeling of the SFR

The results presented in this section can be incorporated into a new subgrid star formation prescription based on the PRFM theory, for implementation in large-scale galaxy formation models where the ISM, star formation, and feedback are unresolved. For this application, S. Hassan et al. (2024, submitted) presents detailed procedures for computing the
<figure><img src="image_16.png" /><figcaption>**Figure 10.** Metallicity and weight dependence of SFR surface density for our full simulation suite. The fit to the TIGRESS-classic suite presented in Ostriker & Kim (2022) is shown as the black dashed line. Symbols and lines have the same meaning as in Figure 8. The observable manifestation of lower feedback efficiency under higher pressure (or density) and higher metallicity conditions is the slightly superlinear dependence of $\Sigma_{SFR}$ on $W$, and the moderate positive dependence of $\Sigma_{SFR}$ on $Z$, quantified in the fit shown.</figcaption></figure>

required quantities from simulation variables, including the calculation of the equilibrium weight $W$ considering the contribution from gaseous and stellar disks as well as a dark matter halo. While the simple weight estimate in Equation (10) is applicable for normal star-forming galaxies where the gas disk is thinner than the stellar disk, the more general weight estimator in S. Hassan et al. (2024, submitted) is applicable for any thickness ratio between gaseous and stellar disks.

Once the weight $W$ is determined, the total feedback yield can be used to obtain the SFR. Within the parameter space explored in this paper, either the single fit to the total feedback yield (Equation (20)) or the sum of fits for thermal (Equation (14)) and turbulent+magnetic (Equation (17)) components can be used to obtain $\Upsilon_{tot}$. While not recommended, if extrapolation beyond the parameter space is unavoidable, using the latter ($\Upsilon_{tot} = \Upsilon_{th} + \Upsilon_{turb+mag}$) is safer as the expected behavior would be captured in low- and high-pressure regimes where thermal and turbulent+magnetic support respectively dominates.

For a disk galaxy, we can express $\Sigma_{SFR}$ in terms of star formation efficiency per dynamical time, $\epsilon_{dyn}$:

$$ \Sigma_{SFR} \equiv \epsilon_{dyn} \frac{\Sigma_{gas}}{t_{dyn}}. \qquad (24) $$

In the above,

$$ t_{dyn} \equiv \frac{2H}{\sigma_{z,eff}} = \frac{2\sigma_{z,eff}}{\langle g_z \rangle} \qquad (25) $$

is the vertical dynamical timescale, where we have used $\Sigma_{gas} = 2\rho_{mid}H$ and $P_{tot} = \rho_{mid}\sigma_{z,eff}^2 = \Sigma_{gas}\langle g_z \rangle/2$. Equation (24) is completely equivalent to Equation (2). Since these expressions relate to gas mass and SFR, they are equally applicable for any mass element that averages spatially over the multiphase ISM and temporally over several cycles of star formation and feedback. We could thus write, for a large-scale, an unresolved mass element $m_{gas}$, a mean SFR

$$ \dot{m}_{\star} = \frac{m_{gas}}{t_{dep}} = \epsilon_{dyn} \frac{m_{gas}}{t_{dyn}}. \qquad (26) $$

Using Equations (11) and (13), the gas depletion time can be expressed in terms of vertical gravity and feedback yield, or alternatively using Equation (25) in terms of the dynamical time:

$$ t_{dep} = \frac{2\Upsilon_{tot}}{\langle g_z \rangle} = \frac{\Upsilon_{tot}}{\sigma_{z,eff}} t_{dyn}. \qquad (27) $$

Section 4.2 shows that the vertical gravity is nearly independent of metallicity, and this can also be seen explicitly from Equations (12) to (23). Because $\langle g_z \rangle$ is essentially independent of metallicity, the first equality in Equation (27) implies that metallicity is expected to affect the timescale for converting gas to stars only through the feedback yield, $\Upsilon_{tot}$. The quantities $\langle g_z \rangle$ and $t_{dyn}$ can be easily estimated from properties of galaxies (in simulations or observations) that are measured on ~kpc scales; Equation (12) gives $\langle g_z \rangle$ for the case where the gas disk is thinner than the stellar disk, and S. Hassan et al. (2024, submitted) present more general algebraic formulae.

From our fits for $\Upsilon_{tot}$ in Equation (20) and $\sigma_{z,eff}$ in Equation (23), we can obtain an expression for the dependence of star formation efficiency on $W$ and $Z'_g$:

$$ \epsilon_{dyn} \equiv \frac{t_{dyn}}{t_{dep}} = \frac{\sigma_{z,eff}}{\Upsilon_{tot}} = 0.0071 W_4^{0.41} Z'_g^{0.30}. \qquad (28) $$

We note that, in this equation, an alternative to using Equation (20) is to use $\Upsilon_{tot} = \Upsilon_{th} + \Upsilon_{turb+mag}$ with the calibrations given in Equations (14) and (17). From the TIGRESS-classic results in Ostriker & Kim (2022), we obtain a similar scaling with $W$ and slightly larger normalization, $\epsilon_{dyn} = 0.012W_4^{0.43}$.

In systems where gas gravity dominates (e.g., self-gravitating clouds), a similar expression to Equation (26) is often used to characterize star formation in terms of an efficiency per gas freefall time (e.g., Krumholz & Tan 2007; Padoan et al. 2014; Evans et al. 2022), with $\epsilon_{dyn} \rightarrow \epsilon_{ff}$ and $t_{dyn} \rightarrow t_{ff} = [3\pi/(32G\bar{\rho}_{gas})]^{1/2}$. The freefall time-based star formation recipe (often called a Schmidt-type star formation recipe; Schmidt 1959) has been widely used in galaxy formation simulations where the multiphase ISM is not resolved (e.g., Somerville & Davé 2015; Naab & Ostriker 2017, for reviews). However, adopting a prescription based solely on the gaseous freefall time on galactic scales, where the gravity from the stellar disk (and sometimes dark matter) is of equal or greater importance, is not well justified.

The advantages of our new subgrid model approach compared to currently used Schmidt-type recipes are thus twofold. (1) The dynamical time estimate (as opposed to gas freefall time) takes into account the gravitational contributions from both gas and external stellar and dark matter gravity. (2) The efficiency per
<figure><img src="image_17.png" /><figcaption>**Figure 11.** Calibration of the effective equation of state; i.e., the relation between total pressure and gas number density averaged over all gas (left) and two-phase gas (right). A slightly shallower relation $P_{\text{tot}} \propto n_{\text{mid}}^{1.29}$ with a weak metallicity dependence is found in the NCR suite compared to that of the TIGRESS-classic ($P_{\text{tot}} \propto n_{\text{mid}}^{1.43}$, black dashed in right panel; Ostriker & Kim 2022). Symbols and lines have the same meaning as in Figure 8.</figcaption></figure>

dynamical time is taken to depend on the galactic environment, via the ISM weight and metallicity, rather than being a fixed constant. In addition, a major advance of this paper (along with Ostriker & Kim 2022) is that the calibrations of $\Upsilon_{\text{tot}}$ and $\sigma_{\text{z,eff}}$ are based on holistic simulations of the ISM with star formation and feedback, with accurate treatments of radiative transfer and photochemistry. This allows us to obtain a physics-based (as opposed to empirical) star formation prescription. Finally, the choice to calibrate $\epsilon_{\text{dyn}}$ as a function of $W$ and $Z'_g$ makes it possible to apply this subgrid prescription even in large-box galaxy formation simulations where the vertical structure of the disk is unresolved. In this situation, $\bar{\rho}_{\text{gas}}$ would be lower than it would be at higher resolution, making $t_{\text{ff}}$ larger than it should be. The underestimation of density at low resolution is due to the limit on density set by the smoothing scale and mass resolution, which may be lower than the true density should be for a given surface density and gravitational potential. However, it is still possible to obtain a robust estimate of $W$ as well as $t_{\text{dyn}}$ from available simulation variables (see S. Hassan et al. 2024, submitted).

# 5. Discussion

## 5.1. When Does Thermal Pressure Become Important?

One of the main results of this paper is the metallicity dependence of SFRs (Figures 3 and 10). We attribute the reduced SFRs at low metallicities to the enhancement of the thermal feedback yield under conditions of reduced FUV attenuation (Section 4.3). Because radiation is much less attenuated in low surface density conditions, $\Upsilon_{\text{th}}$ also increases rapidly with decreasing weight (which depends nonlinearly on surface density). The bivariate fit in Equation (14) is

$\Upsilon_{\text{th}} \propto W^{-0.46} Z_g^{-0.53}$. This means that, at lower metallicities and gas surface densities, we expect more efficient thermal pressure regulation by UV radiative feedback. The turbulent feedback yield $\Upsilon_{\text{turb}}$ depends on the specific momentum injected by SNe (or mechanical feedback in general), which decreases as cooling increases. Since SNRs cool when the shock velocity is slightly higher under conditions of higher ambient density and metallicity (Kim & Ostriker 2015a; Kim et al. 2023b), the turbulent yield follows similar trends to the thermal yield, but with shallower dependencies; combining turbulent with magnetic yields, we find $\Upsilon_{\text{turb+mag}} \propto W^{-0.22} Z_g^{-0.18}$. Because of the different yield scalings, thermal pressure would begin to dominate over nonthermal pressure roughly when $WZ'_g$ falls below $100 \text{ cm}^{-3} \text{ K}$; these combined low-pressure and low-metallicity conditions may be found in dwarfs or the outer parts of spirals.

For the total yield, Equation (20) gives a dependence on weight and metallicity $\Upsilon_{\text{tot}} \propto W^{-0.29} Z_g^{-0.27}$, which compromises between the steeper thermal and shallower nonthermal dependencies. From Equation (27), the increase in total feedback yield at low metallicity would imply a corresponding increase in the gas depletion time $t_{\text{dep}} \propto \Upsilon_{\text{tot}}$, for a given galactic gravitational field. Figure 10 shows our results for the dependence of $\Upsilon_{\text{tot}}$ on $W$ and $Z$ for the full simulation suite. Overall, the enhancement in total feedback yield at $Z' = 0.1$ compared to $Z' = 1$ is about a factor 2–2.5. Given the relatively small predicted variation in $\Sigma_{\text{SFR}}$ or $t_{\text{dep}}$ for varying abundances, and potentially large uncertainties in observational measurements, a quantitative confirmation of these predictions may require large ensemble averages of regions sharing similar conditions.
<figure><img src="image_18.png" /><figcaption>**Figure 12.** Mass-weighted mean vertical velocity dispersions of the 2P phase as defined by Equation (8) for the total (left) and turbulent (right) pressure. On the left, the fit to the TIGRESS-classic suite presented in Ostriker & Kim (2022) is shown as the black dashed line. Symbols and lines have the same meaning as in Figure 8.</figcaption></figure>

While comprehensive observational comparisons are needed, there is already some observational evidence aligned with the predictions of our metallicity-dependent models. Herrera-Camus et al. (2017) have measured thermal pressure in the CNM for nearby galaxies using Herchel [C II] 158 µm observations combined with H I and CO data. The KINGFISH sample they analyzed provides a large data set for >500 atomic-dominated regions with typical sizes of ~1 kpc², showing a statistically significant correlation between thermal pressure and Σ<sub>SFR</sub> (see their Figure 5). Using oxygen abundance measurements from Moustakas et al. (2010), Herrera-Camus et al. (2017) reported a systematic trend of increasing thermal pressure at lower metallicities for a given Σ<sub>SFR</sub>. This is consistent with higher thermal feedback yield at lower metallicities at a given weight predicted in our results (Equation (14)).

Our results have important implications for galaxy formation modeling, where it is difficult to follow the thermal structure of the ISM via explicit treatments of cooling and heating processes. Given the large dynamic range that must be covered in global or cosmological galaxy simulations, it has been challenging to properly capture the thermal state of the ISM. Since turbulent pressure dominates over thermal pressure in the Milky Way and other nearby galaxies, an implicit assumption has often been that an exact recovery of the thermal pressure is not required for accurately predicting the regulation of star formation, with more emphasis instead placed on the momentum injection from SNe that drives turbulent pressure (e.g., Hopkins et al. 2023). However, because thermal pressure is expected to exceed turbulent pressure in low-metallicity, low-pressure environments, care must be taken to ensure that

simulations of these environments include the necessary physical ingredients for realistic heating and cooling.

A subtlety that is often overlooked is that PE heating depends on the electron abundance through the grain charging parameter. It is often assumed in galactic and cosmological simulations that the only source of ionization is UV radiation (either attenuated metagalactic or locally produced UV). In the neutral ISM, FUV is able to ionize weakly bound electrons to make $C^+$, but the primary source of ionization is believed to be H-ionization by low-energy CRs. When CR ionization is not included, the electron fraction is unphysically low for the neutral ISM, leading to PE heating rates that are too low (see Kim et al. 2023b, Section 8.2). This results in a much higher CNM fraction and lower thermal pressure than is realistic. In some simulations, mechanical feedback is boosted by multiplying the terminal radial momentum per SN by a large factor (typically by a factor of 5; Agertz et al. 2013; Semenov et al. 2017; Li et al. 2018). It is possible that this boosting is required to make up for missing thermal support (as well as magnetic support in purely hydrodynamical simulations). While this enhancement of turbulent pressure may produce reasonable SFRs, a concern is that there may be other unintended consequences (e.g., for outflows).

## 5.2. When Do SFRs Depend on Metallicity?

In the set of simulations presented here, with $Z'_g \ge 0.1$, the dominant heating source in the 2P phase is from the PE effect, proportional to dust abundance. At sufficiently low dust abundance, however, the PE heating will fall below CR heating as the dominant term (see, e.g., Figure 6 in both of Bialy & Sternberg 2019; Kim et al. 2023b), and at sufficiently
low dust and metal abundance, the primary cooling in both warm and hot diffuse gas will be from hydrogen (see, e.g., Figures 3, 4, 13; and 1, 3 in Bialy & Sternberg 2019; Kim et al. 2023b, respectively). Thus, sensitivity to metallicity (for both thermal pressure and SN momentum injection) must eventually drop when abundances are low enough. However, the exact point when this occurs depends on the details of the cooling and heating processes, some of which are still uncertain.

Two major pieces of modeling uncertainty affect quantitative predictions for the dominant heating mechanism. First, the dust model adopted in the TIGRESS-NCR framework does not include any metallicity dependence. In reality, metallicity-dependent changes in dust properties would be important for both radiation attenuation and modeling PE heating, which is dominated by small grains and PAHs. There is empirical evidence of metallicity dependence of the PAH abundance, dust-to-metal ratio, and grain size distribution (e.g., Draine et al. 2007; Sandstrom et al. 2012; Rémy-Ruyer et al. 2014; De Vis et al. 2019; Roman-Duval et al. 2022; Chastenet et al. 2023). All of these can cause metallicity dependence of attenuation of UV radiation and the PE heating efficiency (Bakes & Tielens 1994; Weingartner & Draine 2001a). It has been proposed, for example, that the reduction of the PAH abundance at low metallicity may not be gradual but a step-like function at $Z' \lesssim 0.2$ (Draine et al. 2007; but see Aniano et al. 2020). In that case, the overall metallicity dependence we find at $Z' > 0.1$ may still hold, while CR heating would begin to dominate earlier at $Z' \sim 0.2$.

Second, our current treatment of CR ionization rate is still provisional even though it is motivated by empirical relations and first-order physical arguments. The CR transport problem in a multiphase, highly dynamic ISM is not fully understood. Recent advances have been made by using two-moment evolution equations for CR transport in order to capture advection, diffusion, and (at high scattering rate) streaming limited by the Alfvén speed (e.g., Jiang & Oh 2018; Girichidis et al. 2020; Hopkins et al. 2022a). In most studies, the CR scattering rates are imposed as a function of energy, but some implementations (Armillotta et al. 2021; Hopkins et al. 2022b) have considered the self-confinement paradigm (e.g., Zwei- bel 2013) that is believed to hold at $E \lesssim 100$ GeV, in which scattering rates are set by a balance between streaming-driven wave excitation and damping. Provided that the full multiphase structure and dynamics of the ISM are sufficiently resolved by the underlying MHD model (including a hot, high-velocity component that rapidly advects CRs out of the disk, a primarily neutral component where CRs are highly diffusive, and a warm ionized component where transport is limited by Alfvènic streaming), realistic properties of GeV CRs are obtained (Armillotta et al. 2024). An extension of this work to allow for multiple CR energy groups, coupled with the TIGRESS-NCR framework, will be able to shed light on the environmental dependence of both ISM heating from MeV CRs and ISM dynamics driven by GeV CRs.

In addition to uncertainties in the ISM physics, there are also uncertainties related to the initial mass function (IMF) and population synthesis model, especially at low metallicities and high redshifts (e.g., Conroy 2013; Eldridge & Stanway 2022). In this paper, we simply keep both stellar models unchanged from the standard choice we made for solar neighborhood conditions, i.e., STARBURST99 (Leitherer et al. 1999, 2014) coupled with a Kroupa IMF (Kroupa 2001) and the Geneva

evolutionary tracks for nonrotating stars. More realistically, at low metallicities, UV radiation may be stronger for the same population of stars (e.g., Grasha et al. 2021), shifting up the normalization for the thermal feedback yield related to the ratio of FUV radiation to SFR, $\Sigma_{\text{FUV}}/\Sigma_{\text{SFR}}$. Stellar rotation and binary evolution can also change the UV photon production rate as well as SN rate (e.g., Levesque et al. 2012; Eldridge & Stanway 2022). There is also evidence that the IMF may have been more bottom heavy in environments that produce very high stellar densities, based on observations of massive elliptical galaxies (e.g., Conroy & van Dokkum 2012). Further numerical experiments that vary the model ingredients using TIGRESS-NCR or similar frameworks will be critical to understanding the parameter sensitivity in observational diagnostics and underlying physical properties, such as SFRs and ISM properties, with the ultimate goal of constraining the uncertain model parameters.

### 5.3. Comparison with Previous Numerical Work

The most similar numerical simulations to those analyzed here were presented in Hu et al. (2021). They ran a set of simulations representing a kiloparsec patch of solar neighborhood, similar to our $R^8$ series, using the Lagrangian code GIZMO (Hopkins 2015), and varying metallicities over the range $Z' = 0.1-3$. Their simulations include time-dependent hydrogen chemistry, as well as resolved SN feedback similar to TIGRESS-NCR. Taking advantage of high resolution in the dense gas obtained by the Lagrangian approach, they focus on the effect of metallicity in the atomic-to-molecular transition and the distribution of carbon-bearing species that are key observables (C II, C I, and CO).

The long-term mean SFR from the solar metallicity simulation of Hu et al. (2021) is similar to that in our $Z' = 1$ TIGRESS-NCR model, which is also in agreement with observations and with our previous TIGRESS simulations with solar neighborhood conditions (Kim & Ostriker 2017). Although SFRs were not the main focus of their work, it is worth noting that, unlike us, Hu et al. (2021) did not find systematic metallicity dependence of mean SFRs. Since they focused on chemical properties, it is difficult to pinpoint the reason. While the gas surface density and gravitational potential were similar to our $R^8$ models and their heating and cooling prescriptions were similar, a potential key difference is in the treatment of the UV radiation field. Instead of the direct UV radiation transfer approach of this paper, at any time, they first compute a uniform background UV radiation field by scaling proportional to $\Sigma_{\text{SFR}}$, using the observed solar neighborhood radiation field and $\Sigma_{\text{SFR}}$ to normalize. They then apply local dust shielding based on column density, calculated using a tree method with a fixed shielding length of 100 pc. Without direct radiative transfer on large scales, this approach did not capture the overall enhancement in the ratio of the FUV radiation field to the SFR at low metallicity. In our simulations, the reduced large-scale attenuation of UV at low metallicity is what leads to higher thermal feedback yield and ultimately enables a lower SFR.

The difference in outcomes suggests that accounting for large-scale UV radiation transfer may be important when aiming to quantitatively capture self-regulation of SFRs, especially for the regime where thermal pressure plays a major role. If an expensive radiation transfer calculation is not feasible, an effective model calibrated based on more direct
radiation transfer results like ours is an alternative approach. As pointed out in Kim et al. (2023a), the simple plane-parallel approximation proposed in Ostriker et al. (2010; see also Ostriker & Kim 2022) provides a reasonable model for large-scale attenuation of FUV radiation (see also Bialy 2020). A more extensive analysis of radiation fields obtained from ray-tracing in the TIGRESS-NCR suite, and comparison with other approximate models, will be presented in N. Linzer et al. (2024, submitted).

Although their mean SFRs did not depend on metallicity, Hu et al. (2021) observed an increase of burstiness in SFRs as metallicity decreases. We find a similar qualitative trend for both R8 and LGR4 (see Figure 3). The reason for the increase in SFR burstiness in our simulations is twofold: (1) faster quenching of star formation (due to heating of the cold phase) over a larger region at lower Z as the radiation is less attenuated, and (2) slower recovery of cold gas due to longer cooling times from reduced metal cooling.

## 5.4. Caveats and Future Perspectives

We find that the mean metallicity increases by a factor of 1.5–3 over the course of our simulations. However, we ignore this metal enrichment as we do not follow dust evolution explicitly, and we also do not include other aspects of global evolution (such as changing gas surface density from accretion). Hu et al. (2023b) included a model of dust formation and destruction in their simulations and found that dust growth is fast and reaches high dust abundance quickly in dense, star-forming gas. We cover a range of dust abundance at low metallicity, but we did not consider spatially varying dust abundance. In the future, it would be interesting to couple our TIGRESS-NCR framework to dust evolution models and investigate the effect of locally varying metallicity, dust abundance, and dust-to-metal ratios.

Finally, it is worth mentioning that the regime of high weight/pressure ($W/k_B \gtrsim 10^{5.5} \text{ cm}^{-3} \text{ K}$) is still subject to some uncertainty, as evident in the results presented here as well as other ongoing studies. In particular, the saturation level of magnetic fields and turbulence properties in the high-weight/-pressure regime may have some sensitivity to the simulation box size as it sets the largest scale of collapse and hence clustering of SNe. In addition, at low redshift, high-weight/-pressure regimes are usually found within distinct galactic structures such as spiral arms (Kim et al. 2020c) and nuclear rings fed by bars (Moon et al. 2023). To reach firm conclusions regarding the quantitative results for turbulent and magnetic feedback yields in the high-weight/-pressure regime and the effective equation of state, further investigation is warranted. Recent efforts in coupling various photochemistry treatments with radiation transfer (usually based on the two-moment method with the M1 closure; e.g., Kannan et al. 2020b; Lupi et al. 2020; Katz et al. 2022) in cosmological and global isolated galaxy modeling will help shed light on exploring more extreme conditions. Nevertheless, box size and global geometry would not alter the robustness of our conclusions regarding the qualitative behavior of feedback yields with metallicity and weight, or the insensitivity of vertical dynamical equilibrium and effective equation of state to metallicity.

# 6. Summary

In this paper, we present the first results from a new suite of local galactic patch simulations with varying galactic conditions (weight and metallicity) using the TIGRESS-NCR framework (Kim et al. 2023a). The new simulation suite includes a total of 28 models covering a wide range of gas and stellar surface densities $\Sigma_{\text{gas}} = 5-150 \text{ M}_\odot \text{ pc}^{-2}$ and $\Sigma_{\text{star}} = 1-50 \text{ M}_\odot \text{ pc}^{-2}$. We vary gas metallicity and dust abundance up to 3 times solar values, metals down to 0.1 times solar, and dust down to 0.025 times solar. After running for at least one orbit time to pass through an early transient stage, we run all models for more than two orbit times, corresponding to evolution periods ranging from 150 Myr to 1.5 Gyr (shorter duration for higher gas and SFR surface density conditions). The emergent midplane total pressure from our suite covers a range $P_{\text{tot}}/k_B = 10^3-10^6 \text{ cm}^{-3} \text{ K}$, while the range of SFR surface density is $\Sigma_{\text{SFR}} = 10^{-4}-0.5 \text{ M}_\odot \text{ kpc}^{-2} \text{ yr}^{-1}$.

The TIGRESS-NCR framework represents a significant advance from the original TIGRESS framework (Kim & Ostriker 2017, referred to as “TIGRESS-classic”). The key improvements over TIGRESS-classic are explicit UV radiation transfer using ART (Kim et al. 2017b) from star cluster particles for both nonionizing and ionizing UV, and a photochemistry module, which enables realistic ISM cooling and heating over a range of metal and dust abundances (Kim et al. 2023b). TIGRESS-classic adopted much simpler temperature-dependent cooling and spatially constant (but SFR-dependent) heating, and was only applicable for solar metallicity. Other major dynamical and feedback processes—including MHD, galactic sheared rotation, and self-gravity, as well as star cluster formation in gravitational collapse, and resolved SN feedback—are the same in TIGRESS-classic and TIGRESS-NCR.

The primary focus of this paper is to investigate the metallicity dependence of SFRs in the context of the PRFM star formation theory (e.g., Ostriker & Kim 2022). A key motivation is to calibrate subgrid star formation models that may be implemented in cosmological galaxy formation simulations to follow evolution over a range of redshifts.

Our main findings are as follows:

1. When other conditions are fixed, SFRs drop at lower metallicities (Section 3.2). From $Z' = 1$ to $Z' = 0.1$, we find that the mean $\Sigma_{\text{SFR}}$ is reduced by a factor of 2–3. In terms of gas depletion time $t_{\text{dep}} = \Sigma_{\text{gas}}/\Sigma_{\text{SFR}}$, we find a comparable or slightly larger enhancement in mean $t_{\text{dep}}$ at lower metallicities. This is because reduced $\Sigma_{\text{SFR}}$ in the early evolution results in slightly higher $\Sigma_{\text{gas}}$ in the later evolution. At $Z'_g = 0.1$, a reduction in dust abundance from $Z'_d = 0.1$ to $Z'_d = 0.025$ does not cause a further reduction in $\Sigma_{\text{SFR}}$.

2. In all simulations, the vertical dynamical equilibrium is satisfied (Section 4.2). This means that, at the disk midplane, the total pressure $P_{\text{tot}}$ matches the weight per unit area $W$; vertical profiles of pressure, and weights also match each other within $|z| \lesssim 500 \text{ pc}$. The simple dynamical equilibrium estimator (Equation (10)) is in good agreement with the actual weight irrespective of metallicity. Since the weight is determined mainly by $\Sigma_{\text{gas}}$ and $\Sigma_{\text{star}}$ and is insensitive to metallicity, the total midplane pressure is also insensitive to metallicity. However, the relative contributions to $P_{\text{tot}}$ of the thermal
pressure ($P_{\text{th}}$), the vertical component of Reynolds stress ($P_{\text{turb}} = \rho v_z^2$), and the vertical component of Maxwell stress ($\Pi_{\text{mag}} = B^2/8\pi - B_z^2/4\pi$) vary with metallicity. The increase of $P_{\text{th}}/P_{\text{tot}}$ at low metallicity (by more than a factor of 2 from $Z'_g = 3$ to $Z'_g = 0.1$ in the R8 series) implies a corresponding reduction in the other fractional pressure contributions. The hot, WIM, and 2p phases are in total pressure equilibrium within the scale height of the 2p phase.

3. We provide a bi-variate power-law fit for each feedback yield component and the total feedback yield, as functions of weight and metallicity (Section 4.3). We calculate the mean values of the stresses and SFR surface density over an extended interval during the later evolution covering multiple star formation and feedback cycles. We then measure the ratios of each stress component to the SFR surface density, $\Sigma_{\text{SFR}}$. These ratios are termed feedback yields in the PRFM theory, denoted $\Upsilon_{\text{tot}} = P_{\text{tot}}/\Sigma_{\text{SFR}}$, $\Upsilon_{\text{turb}} = P_{\text{turb}}/\Sigma_{\text{SFR}}$, $\Upsilon_{\text{th}} = P_{\text{th}}/\Sigma_{\text{SFR}}$, and $\Upsilon_{\text{mag}} = \Pi_{\text{mag}}/\Sigma_{\text{SFR}}$. We confirm, consistent with Ostriker & Kim (2022), a decreasing trend in all feedback yields as the weight increases. We find a clear increasing trend in $\Upsilon_{\text{th}}$ with decreasing metallicity down to $Z'_g = 0.1$. There is also a strong increase in $\Upsilon_{\text{th}}$ with decreasing weight because environments with low pressure ($P_{\text{tot}} \approx W$) have low density such that radiation is only weakly attenuated and efficiently heats the gas. There is almost no metallicity dependence in $\Upsilon_{\text{turb}}$. While the individual behavior of turbulent and magnetic feedback yields becomes irregular at high weights ($W/k_B > 10^{5.5} \text{ cm}^{-3} \text{ K}$), the combined nonthermal (turbulent+magnetic) feedback yield $\Upsilon_{\text{turb+mag}}$ is better described by a power law with overall reduced scatter. $\Upsilon_{\text{turb+mag}}$ depends less sensitively on both weight and metallicity, compared to $\Upsilon_{\text{th}}$.

4. We provide a new calibration for the effective equation of state for star-forming multiphase gas (Section 4.4). We obtain a shallower exponent than TIGRESS-classic (1.29 versus 1.43 from Ostriker & Kim 2022) for the total pressure and density relation with no metallicity dependence. We also provide a fit for the effective vertical velocity dispersion of the warm-cold 2p medium, which increases weakly at higher pressure but is insensitive to metallicity.

5. We describe how our calibrations of feedback yield and effective equation of state can be used in a new subgrid star formation model based on the PRFM theory (Section 4.5). This model expresses the SFR in terms of the dynamical time $t_{\text{dyn}}$ and an efficiency factor $\epsilon_{\text{dyn}}$ (Equation (26)). We provide a fit for $\epsilon_{\text{dyn}}$ based on our calibrations for $\Upsilon_{\text{tot}}$ and $\sigma_{z,\text{eff}}$ in terms of weight $W$ and metallicity (Equation (28)). The advantages of this new model are that it is grounded in explicit simulations of the ISM with accurate physics treatments and widely ranging galactic environments and that it takes into account the gravity of stars and dark matter as well as gas. By using weight rather than density as the input parameter, our SFR model is designed to be robustly applicable in large-box galaxy formation simulations in which gas disk scale heights are unresolved (and hence, mean density is resolution dependent).

The main focus of this work was on studying the variation of SFRs under a wide range of galactic conditions, covering

environments observed in both massive galaxies and dwarfs, and in the low- and high-redshift universe. Here, we employed the “feedback yield” to quantitatively characterize the equilibrium multiphase ISM pressure response to energy inputs from recently formed stars. An equally interesting numerical problem, and equally important for developing new cosmological subgrid models, is to quantify how galactic multiphase outflows are produced due to star formation feedback. Following Kim et al. (2020b, 2020a), galactic outflows can be characterized in terms of overall loading factors as well as distributions of outflow speed and sound speed. The present simulation suite, and additional TIGRESS-NCR simulations with augmented physics and parameter coverage, will enable extensions of previous wind-driving analyses to cover a range of metallicities.

## Acknowledgments

We acknowledge the anonymous referee for a useful report. We are grateful to Amiel Sternberg for constructive comments on the manuscript. This work was supported by grant 10013948 from the Simons Foundation to Princeton University, sponsoring the Learning the Universe collaboration. The work of C.-G.K. was partly supported by NASA ATP grant No. 80NSSC22K0717. J.-G.K. acknowledges support from the EACOA Fellowship awarded by the East Asia Core Observatories Association. G.L.B. acknowledges support from the NSF (AST-2108470 and AST-2307419, ACCESS), a NASA TCAN award, and the Simons Foundation. S.H. acknowledges support for program No. HST-HF2-51507 provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, incorporated, under NASA contract NAS5-26555. The Center for Computational Astrophysics at the Flatiron Institute is supported by the Simons Foundation. Resources supporting this work were provided in part by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and in part by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center. This research has made use of NASA’s Astrophysics Data System.

*Software:* Athena (Stone et al. 2008; Stone & Gardner 2009), astropy (Astropy Collaboration et al. 2013, 2018, 2022), scipy (Virtanen et al. 2020), numpy (van der Walt et al. 2011), IPython (Perez & Granger 2007), matplotlib (Hunter 2007), xarray (Hoyer & Hamman 2017), pandas (McKinney 2010), CMasser (van der Velden 2020), adstex (https://github.com/yymao/adstex).

## ORCID iDs

Chang-Goo Kim https://orcid.org/0000-0003-2896-3725

Eve C. Ostriker https://orcid.org/0000-0002-0509-9113

Jeong-Gyu Kim https://orcid.org/0000-0001-6228-8634

Munan Gong https://orcid.org/0000-0003-1613-6263

Greg L. Bryan https://orcid.org/0000-0003-2630-9228

Drummond B. Fielding https://orcid.org/0000-0003-3806-8548

Sultan Hassan https://orcid.org/0000-0002-1050-7572

Matthew Ho https://orcid.org/0000-0003-3207-8868

Sarah M. R. Jeffreson https://orcid.org/0000-0002-4232-0200
Rachel S. Somerville https://orcid.org/0000-0002-6748-6821

Ulrich P. Steinwandel https://orcid.org/0000-0001-8867-5026

References

Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25

Aniano, G., Draine, B. T., Hunt, L. K., et al. 2020, ApJ, 889, 150

Armillotta, L., Ostriker, E. C., & Jiang, Y.-F. 2021, ApJ, 922, 11

Armillotta, L., Ostriker, E. C., Kim, C.-G., & Jiang, Y.-F. 2024, ApJ, 964, 99

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822

Barrera-Ballesteros, J. K., Sánchez, S. F., Heckman, T., et al. 2021, MNRAS, 503, 3643

Bialy, S. 2020, ApJ, 903, 62

Bialy, S., & Sternberg, A. 2019, ApJ, 881, 160

Bieri, R., Naab, T., Geen, S., et al. 2023, MNRAS, 523, 6336

Boulares, A., & Cox, D. P. 1990, ApJ, 365, 544

Buck, T., Obreja, A., Macciò, A. V., et al. 2020, MNRAS, 491, 3461

Chastenet, J., Sutter, J., Sandstrom, K., et al. 2023, ApJL, 944, L11

Chevance, M., Krumholz, M. R., McLeod, A. F., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 1

Colling, C., Hennebelle, P., Geen, S., Iffrig, O., & Bournaud, F. 2018, A&A, 620, A21

Conroy, C. 2013, ARA&A, 51, 393

Conroy, C., & van Dokkum, P. G. 2012, ApJ, 760, 71

Cox, D. P. 1972, ApJ, 178, 159

De Vis, P., Jones, A., Viaene, S., et al. 2019, A&A, 623, A5

Deng, Y., Li, H., Kannan, R., et al. 2024, MNRAS, 527, 478

Draine, B. T. 1978, ApJS, 36, 595

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton, NJ: Princeton Univ. Press)

Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

El-Badry, K., Ostriker, E. C., Kim, C.-G., Quataert, E., & Weisz, D. R. 2019, MNRAS, 490, 1961

Eldridge, J. J., & Stanway, E. R. 2022, ARA&A, 60, 455

Evans, N. J., Kim, J.-G., & Ostriker, E. C. 2022, ApJL, 929, L18

Feldmann, R., Quataert, E., Faucher-Giguère, C.-A., et al. 2023, MNRAS, 522, 3831

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, RMXA, 53, 385

Field, G. B., Goldsmith, D. W., & Habing, H. J. 1969, ApJL, 155, L149

Fielding, D., Quataert, E., & Martizzi, D. 2018, MNRAS, 481, 3325

Gatto, A., Walch, S., Naab, T., et al. 2017, MNRAS, 466, 1903

Geen, S., Bieri, R., Rosdahl, J., & de Koter, A. 2021, MNRAS, 501, 1352

Geen, S., Hennebelle, P., Tremblin, P., & Rosdahl, J. 2016, MNRAS, 463, 3129

Geen, S., Pellegrini, E., Bieri, R., & Klessen, R. 2020, MNRAS, 492, 915

Gentry, E. S., Krumholz, M. R., Madau, P., & Lupi, A. 2019, MNRAS, 483, 3647

Girard, M., Fisher, D. B., Bolatto, A. D., et al. 2021, ApJ, 909, 12

Girichidis, P., Naab, T., Hanasz, M., & Walch, S. 2018, MNRAS, 479, 3042

Girichidis, P., Pfrommer, C., Hanasz, M., & Naab, T. 2020, MNRAS, 491, 993

Girichidis, P., Walch, S., Naab, T., et al. 2016, MNRAS, 456, 3432

Gnat, O., & Ferland, G. J. 2012, ApJS, 199, 20

Gnat, O., & Sternberg, A. 2007, ApJS, 168, 213

Grasha, K., Roy, A., Sutherland, R. S., & Kewley, L. J. 2021, ApJ, 908, 241

Haffner, L. M., Dettmar, R. J., Beckman, J. E., et al. 2009, RvMP, 81, 969

Haid, S., Walch, S., Seifried, D., et al. 2018, MNRAS, 478, 4799

Herrera-Camus, R., Bolatto, A., Wolfire, M., et al. 2017, ApJ, 835, 201

Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555

Hopkins, P. F. 2015, MNRAS, 450, 53

Hopkins, P. F., Butsky, I. S., Panopoulou, G. V., et al. 2022a, MNRAS, 516, 3470

Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581

Hopkins, P. F., Squire, J., Butsky, I. S., & Ji, S. 2022b, MNRAS, 517, 5413

Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800

Hopkins, P. F., Wetzel, A., Wheeler, C., et al. 2023, MNRAS, 519, 3154

Hoyer, S., & Hamman, J. 2017, JOSS, 5, 10

Hu, C.-Y., Naab, T., Glover, S. C. O., Walch, S., & Clark, P. C. 2017, MNRAS, 471, 2151

Hu, C.-Y., Smith, M. C., Teyssier, R., et al. 2023a, ApJ, 950, 132

Hu, C.-Y., Sternberg, A., & van Dishoeck, E. F. 2021, ApJ, 920, 44

Hu, C.-Y., Sternberg, A., & van Dishoeck, E. F. 2023b, ApJ, 952, 140

Hu, C.-Y., Zhukovska, S., Somerville, R. S., & Naab, T. 2019, MNRAS, 487, 3252

Hunter, J. D. 2007, CSE, 9, 90

Iffrig, O., & Hennebelle, P. 2015, A&A, 576, A95

Indriolo, N., Neufeld, D. A., Gerin, M., et al. 2015, ApJ, 800, 40

Jeffreson, S. M. R., Krumholz, M. R., Fujimoto, Y., et al. 2021, MNRAS, 505, 3470

Jeffreson, S. M. R., Sun, J., & Wilson, C. D. 2022, MNRAS, 515, 1663

Jiang, Y.-F., & Oh, S. P. 2018, ApJ, 854, 5

Joung, M. K. R., & Mac Low, M.-M. 2006, ApJ, 653, 1266

Joung, M. R., Mac Low, M.-M., & Bryan, G. L. 2009, ApJ, 704, 137

Kannan, R., Marinacci, F., Simpson, C. M., Glover, S. C. O., & Hernquist, L. 2020a, MNRAS, 491, 2088

Kannan, R., Marinacci, F., Vogelsberger, M., et al. 2020b, MNRAS, 499, 5732

Karpov, P. I., Martizzi, D., Macias, P., et al. 2020, ApJ, 896, 66

Katz, H. 2022, MNRAS, 512, 348

Katz, H., Liu, S., Kimm, T., et al. 2022, arXiv:2211.04626

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Kim, C.-G., Kim, J.-G., Gong, M., & Ostriker, E. C. 2023a, ApJ, 946, 3

Kim, C.-G., Kim, W.-T., & Ostriker, E. C. 2011, ApJ, 743, 25

Kim, C.-G., & Ostriker, E. C. 2015a, ApJ, 802, 99

Kim, C.-G., & Ostriker, E. C. 2015b, ApJ, 815, 67

Kim, C.-G., & Ostriker, E. C. 2017, ApJ, 846, 133

Kim, C.-G., & Ostriker, E. C. 2018, ApJ, 853, 173

Kim, C.-G., Ostriker, E. C., Fielding, D. B., et al. 2020a, ApJL, 903, L34

Kim, C.-G., Ostriker, E. C., & Raileanu, R. 2017a, ApJ, 834, 25

Kim, C.-G., Ostriker, E. C., Somerville, R. S., et al. 2020b, ApJ, 900, 61

Kim, J.-G., Gong, M., Kim, C.-G., & Ostriker, E. C. 2023b, ApJS, 264, 10

Kim, J.-G., Kim, W.-T., Ostriker, E. C., & Skinner, M. A. 2017b, ApJ, 851, 93

Kim, J.-G., Kim, W.-T., & Ostriker, E. C. 2018, ApJ, 859, 68

Kim, J.-G., Ostriker, E. C., & Filippova, N. 2021, ApJ, 911, 128

Kim, W.-T., Kim, C.-G., & Ostriker, E. C. 2020c, ApJ, 898, 35

Kimm, T., & Cen, R. 2014, ApJ, 788, 121

Kroupa, P. 2001, MNRAS, 322, 231

Krumholz, M. R., McKee, C. F., & Bland-Hawthorn, J. 2019, ARA&A, 57, 227

Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304

Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 571

Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021a, ApJ, 914, 90

Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021b, ApJL, 922, L3

Lee, J., Kimm, T., Blaizot, J., et al. 2022, ApJ, 928, 144

Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Levesque, E. M., Leitherer, C., Ekstrom, S., Meynet, G., & Schaerer, D. 2012, ApJ, 751, 67

Li, C., Li, H., Cui, W., et al. 2024, MNRAS, 529, 4073

Li, H., Gnedin, O. Y., & Gnedin, N. Y. 2018, ApJ, 861, 107

Li, H., Vogelsberger, M., Marinacci, F., Sales, L. V., & Torrey, P. 2020, MNRAS, 499, 5862

Lupi, A., Pallottini, A., Ferrara, A., et al. 2020, MNRAS, 496, 5160

Martizzi, D., Faucher-Giguère, C.-A., & Quataert, E. 2015, MNRAS, 450, 504

Martizzi, D., Fielding, D., Faucher-Giguère, C.-A., & Quataert, E. 2016, MNRAS, 459, 2311

McKee, C. F., & Ostriker, J. P. 1977, ApJ, 218, 148

McKinney, W. 2010, in Proc. of the 9th Python in Science Conf., ed. S. van der & J. Millman, 56

Moon, S., Kim, W.-T., Kim, C.-G., & Ostriker, E. C. 2021, ApJ, 914, 9

Moon, S., Kim, W.-T., Kim, C.-G., & Ostriker, E. C. 2023, ApJ, 946, 114

Moustakas, J., Kennicutt, R. C. J., Tremonti, C. A., et al. 2010, ApJS, 190, 233

Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59

Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234

Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163

Oku, Y., Tomida, K., Nagamine, K., Shimizu, I., & Cen, R. 2022, ApJS, 262, 9

Ostriker, E. C., & Kim, C.-G. 2022, ApJ, 936, 137

Ostriker, E. C., McKee, C. F., & Leroy, A. K. 2010, ApJ, 721, 975

Ostriker, E. C., & Shetty, R. 2011, ApJ, 731, 41

Ostriker, J. P., & McKee, C. F. 1988, RvMP, 60, 1
Padoan, P., Federrath, C., Chabrier, G., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. of Arizona Press), 77

Perez, F., & Granger, B. E. 2007, CSE, 9, 21

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Rathjen, T.-E., Naab, T., Girichidis, P., et al. 2021, MNRAS, 504, 1039

Rathjen, T.-E., Naab, T., Walch, S., et al. 2023, MNRAS, 522, 1843

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31

Roman-Duval, J., Jenkins, E. B., Tchernyshyov, K., et al. 2022, ApJ, 928, 90

Sandstrom, K. M., Bolatto, A. D., Bot, C., et al. 2012, ApJ, 744, 20

Schaye, J., & Dalla Vecchia, C. 2008, MNRAS, 383, 1210

Schinnerer, E., Leroy, A., Blanc, G., et al. 2019, Msngr, 177, 36

Schmidt, M. 1959, ApJ, 129, 243

Schneider, E. E., & Mao, S. A. 2024, ApJ, 966, 37

Schneider, E. E., Ostriker, E. C., Robertson, B. E., & Thompson, T. A. 2020, ApJ, 895, 43

Semenov, V. A., Kravtsov, A. V., & Gnedin, N. Y. 2017, ApJ, 845, 133

Smith, M. C., Bryan, G. L., Somerville, R. S., et al. 2021, MNRAS, 506, 3882

Smith, M. C., Fielding, D. B., Bryan, G. L., et al. 2024, MNRAS, 527, 1216

Smith, R. J., Treß, R. G., Sormani, M. C., et al. 2020, MNRAS, 492, 1594

Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289

Steinwandel, U. P., Bryan, G. L., Somerville, R. S., Hayward, C. C., & Burkhart, B. 2023, MNRAS, 526, 1408

Steinwandel, U. P., & Goldberg, J. A. 2023, arXiv:2310.11495

Steinwandel, U. P., Kim, C.-G., Bryan, G. L., et al. 2024, ApJ, 960, 100

Steinwandel, U. P., Moster, B. P., Naab, T., Hu, C.-Y., & Walch, S. 2020, MNRAS, 495, 1035

Stone, J. M., & Gardiner, T. 2009, NewA, 14, 139

Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142

Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137

Sun, J., Leroy, A. K., Ostriker, E. C., et al. 2020, ApJ, 892, 148

Sun, J., Leroy, A. K., Ostriker, E. C., et al. 2023, ApJL, 945, L19

Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253

Tan, B., & Fielding, D. B. 2024, MNRAS, 527, 9683

Thornton, K., Gaudlitz, M., Janka, H. T., & Steinmetz, M. 1998, ApJ, 500, 95

Treß, R. G., Sormani, M. C., Smith, R. J., et al. 2021, MNRAS, 505, 5438

van der Velden, E. 2020, JOSS, 5, 204

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Vandenbroucke, B., & Wood, K. 2018, A&C, 23, 40

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, NatRP, 2, 42

Walch, S., Girichidis, P., Naab, T., et al. 2015, MNRAS, 454, 238

Walch, S., & Naab, T. 2015, MNRAS, 451, 2757

Watson, W. D. 1972, ApJ, 176, 103

Weingartner, J. C., & Draine, B. T. 2001a, ApJS, 134, 263

Weingartner, J. C., & Draine, B. T. 2001b, ApJ, 548, 296

Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99

Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152

Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278

Wolfire, M. G., Vallini, L., & Chevance, M. 2022, ARA&A, 60, 247

Zweibel, E. G. 2013, PhPl, 20, 055501