Introducing TIGRESS-NCR: I. Co-Regulation of the Multiphase Interstellar Medium and Star Formation Rates

Massive, young stars are the main source of energy that maintains multiphase structure and turbulence in the interstellar medium (ISM), and without this"feedback"the star formation rate (SFR) would be much higher than is observed. Rapid energy loss in the ISM and efficient energy recovery by stellar feedback lead to co-regulation of SFRs and the ISM state. Realistic approaches to this problem should solve the dynamical evolution of the ISM, including star formation, and the input of feedback energy self-consistently and accurately. Here, we present the TIGRESS-NCR numerical framework, in which UV radiation, supernovae, cooling and heating processes, and gravitational collapse are modeled explicitly. We use an adaptive ray tracing method for UV radiation transfer from star clusters represented by sink particles, accounting for attenuation by dust and gas. We solve photon-driven chemical equations to determine the abundances of H (time-dependent) and C/O-bearing species (steady-state), which then set cooling and heating rates self-consistently. Applying these methods, we present high-resolution magnetohydrodynamics simulations of differentially rotating local galactic disks representing typical conditions of nearby star-forming galaxies. We analyze ISM properties and phase distributions and show good agreement with existing multiwavelength galactic observations. We measure midplane pressure components (turbulent, thermal, and magnetic) and the weight, demonstrating that vertical dynamical equilibrium holds. We quantify the ratios of pressure components to the SFR surface density, which we call the feedback yields. The TIGRESS-NCR framework will allow for a wide range of parameter exploration, including low metallicity system.


INTRODUCTION
The interstellar medium (ISM) is at the core of the ecosystem in star-forming galaxies. The ISM gives birth to stars and also processes the energy and metals these stars produce, using the majority in maintaining the ISM state while expelling a fraction to larger scales. Modeling the ISM is fundamental to astrophysics, but challenging for many reasons.
Proper treatment of ISM dynamics and energetics must involve simultaneous modeling of the formation and evolution of massive young stars, encompassing the physics that controls star formation, i.e., gravity, tur-cgkim@astro.princeton.edu bulence, and magnetic fields (McKee & Ostriker 2007). The ISM is highly dissipative, quickly losing energy through radiative processes (e.g., Spitzer 1978;Draine 2011). Without continuous and efficient energy inputs, rapid gravitational collapse (approaching the free fall rate) would occur, which would produce far higher star formation rates (SFRs) than those observed (e.g., Sun et al. 2022). Stars can provide the necessary energy, with UV radiation and supernovae (SNe) from massive young stars the two most energetically dominant channels (e.g., Leitherer et al. 1999).
UV radiation is the primary driver of key thermal and chemical processes in warm and cold ISM phases, setting cooling and heating rates (Wolfire et al. 2022). The gas thermal pressure in warm and cold ISM phases depends on the balance of heating and cooling (Field et al. 1969;Wolfire et al. 1995), and at a given temperature, the chemical state of the most abundant atom, hydrogen, can vary significantly: the warm medium (T ∼ 10 4 K) can be neutral or ionized, and the cold medium (T ∼ 10 2 K) can be neutral or molecular. Thermal and chemical states depend on the strength of the UV radiation field as well as the cosmic ray (CR) rate (Wolfire et al. 2003), and these vary spatially due to proximity of sources and shielding by the highly inhomogeneous structure (e.g., Peters et al. 2017). Additionally, shocks driven by SNe both accelerate and heat the gas (Cox 1972;Cioffi et al. 1988). Given the galactic SN rate, the hot medium (T ∼ 10 6 K) created by SN shocks can occupy a large volume (Cox & Smith 1974;McKee & Ostriker 1977) and break out the disk (Shapiro & Field 1976). Turbulence in the warm and cold ISM is driven by the interaction with expanding supernovaheated bubbles (e.g., Korpi et al. 1999;Joung et al. 2009;Kim et al. 2013;Hennebelle & Iffrig 2014).
From many decades of research, a consensus has been reached that multiple distinct phases coexist in the ISM, spanning a wide range of temperature and density (e.g., Ferrière 2001;Cox 2005). Furthermore, because the ISM is dynamic, the thermal and chemical states of any given gas parcel are transient, and states that would traditionally be considered unstable are continuously repopulated. Extensive observational surveys of the multiphase ISM have detailed the properties of the cold molecular gas (Dame et al. 2001;Heyer & Dame 2015); the cold and warm atomic gas (HI4PI Collaboration et al. 2016;Peek et al. 2018), with evidence of a significant fraction being in the unstable temperature regime (Heiles & Troland 2003;Murray et al. 2018); the warm ionized medium (Hill et al. 2008;Krishnarao et al. 2017); and the hot ionized medium (Cox & Reynolds 1987;Bowen et al. 2008).
The existence of the multiphase ISM and the ubiquitous high-level turbulence in it are clear evidence that stellar feedback energy is effectively coupled to the ISM. Feedback leads to inefficient star formation in terms of gas consumption, resulting in observed SFRs that are two orders of magnitude lower than the free-fall rate (Utomo et al. 2018). Because energy dissipation leads to localized collapse and star formation, while the bulk ISM's energy loss can be efficiently recovered from stel-lar feedback, 1 the ISM's physical state and the SFR are intimately connected and in fact must be co-regulated (Ostriker et al. 2010).
To quantitatively understand the co-regulation process in the star-forming ISM, a holistic approach is required, using direct numerical simulations to solve the magnetohydrodynamics (MHD) equations with gravity and explicit cooling and heating. These numerical simulations explicitly model the ISM in all phases, while tracking star formation in gravitationally collapsing regions and directly following energy inputs from recently formed stars. In order to be self-consistent, all gas heating and turbulence driving must either originate in energy inputs from stars, or develop as a consequence of naturally-occurring ISM dynamics (driven by gravity, shear, etc). Realistic numerical simulations of the starforming multiphase ISM require both high mass and spatial resolution to resolve both the mass-dominating (cold and warm) and volume-filling (warm and hot) components. Given the strict resolution requirements for multiphase ISM simulations, to date, the outer dimensions of simulation domains are still limited to a few kpc, corresponding either to low-mass dwarf galaxies (Hu et al. 2017;Steinwandel et al. 2022) or to portions of larger galactic disks, using vertically-stratified boxes (Kim & Ostriker 2017;Gatto et al. 2017;Peters et al. 2017;Rathjen et al. 2021;Kannan et al. 2020).
The TIGRESS framework was developed by Kim & Ostriker (2017) to study the star-forming multiphase ISM, including a full complement of physics modules, sufficient resolution to follow key processes, and computational performance that enables both long-term evolution and comparative study of multiple galactic environments. In the TIGRESS framework, the MHD equations in a local patch of a differentially rotating galactic disk are solved with the grid-based code Athena Stone & Gardiner 2009). Self-gravity of gas is included, together with a fixed vertical potential from stars and dark matter. Cooling is modeled by a temperature-dependent tabulated cooling function appropriate for solar metallicity (Sutherland & Dopita 1993;Koyama & Inutsuka 2002). Sink particles representing star clusters are introduced within cells undergoing runaway gravitational collapse (Gong & Ostriker 2013). Photoelectric heating by FUV radiation is set to scale with the globally attenuated FUV luminosity from star clusters formed in the simulation. Explosions of individual SNe are directly modeled, resolving the Sedov-Taylor stage during which the radial momentum of expanding bubbles is built up and the hot ISM is created in strong shocks (Kim & Ostriker 2015a). Systems modeled by the TIGRESS framework successfully evolve to a quasi-steady state over many star formation and feedback cycles. A large suite of individual TIGRESS simulations covering varying galactic conditions shows that in all cases a state with self-consistently regulated SFR and ISM state is reached (Ostriker & Kim 2022), with multiphase outflows launched from the disk Vijayan et al. 2020;Kim et al. 2020a). These TIGRESS simulations have also been used to study cloud and star formation correlations (Mao et al. 2020), molecular chemistry (Gong et al. 2018(Gong et al. , 2020, diffuse ionized gas (Kado-Fong et al. 2020), polarized dust emission , and CR transport (Armillotta et al. 2021(Armillotta et al. , 2022. In addition, the TIGRESS computational framework of Kim & Ostriker (2017) has been extended to simulate regions with strong spiral structure (Kim et al. 2020b), nuclear rings where bar-driven gas inflows accumulate (Moon et al. 2021(Moon et al. , 2022, and ram pressure stripping by the intracluster medium (Choi et al. 2022).
This paper presents the first application of an extension of the TIGRESS framework, called TIGRESS-NCR, where "NCR" stands for Non-equilibrium Cooling and Radiation.
The two salient new features of TIGRESS-NCR are explicit UV radiation transfer using the adaptive ray tracing method implemented in Athena (Kim et al. 2017b) and the photochemistry model introduced by Kim et al. (2023). We solve time-dependent chemical rate equations for hydrogen species and obtain other abundances assuming formation-destruction balance given the hydrogen species abundances. Cooling in the cold and warm medium in TIGRESS-NCR is determined by abundances of hydrogen species (H, H + , H 2 ) and major coolants (C + , C, CO, O, O + ). Cooling in warm ionized gas is treated with a nebular cooling function that assumes a fixed abundance pattern characteristic of photoionized regions. High-temperature He and metal cooling assume collisional ionization equilibrium (CIE). To follow UV radiation, photon packets emanating from star clusters are transferred along rays through the ISM, with absorption by dust and gas. The major radiative heating channels (photoelectric and photoionization heating) and expansion of overpressurized H II regions (driven by photoionized gas and radiation pressure) are directly modeled. We also include cosmic-ray induced ionization and heating with an attenuation factor inversely proportional to an effective mean column density. Our new framework with adaptive ray tracing improves upon the accuracy of radiation transfer solutions compared to other ISM simulations that adopt more approximate methods, including the local attenuation and local ionization approach in Hu et al. (2021), the tree-based backward radiation transfer method in Rathjen et al. (2021), and the two-moment approach with M1 closure in Kannan et al. (2020); Katz et al. (2022).
In this paper, we focus on technical aspects of the TIGRESS-NCR implementation, and demonstrate how the ISM state and SFR are co-regulated by the full physics treatments in the TIGRESS-NCR framework. We consider two different galactic conditions for the models in this paper, one similar to the solar neighborhood, and one representing inner-galaxy environments. In a companion paper, we use the TIGRESS-NCR implementation to conduct a set of controlled numerical experiments in which we turn on and off individual feedback channels and the magnetic field, in order to investigate the role of each physical process in regulating SFRs and ISM properties. In subsequent papers, we will present detailed analyses of radiation properties, ISM energetics, and galactic outflows.
We describe numerical details in Section 2, drawing from Kim & Ostriker (2017) and from Kim et al. (2023). TIGRESS-NCR specific treatments regarding truncation of rays for efficient calculations are detailed in Section 2.3 and Appendix A. Section 3 describes the ISM properties, energetics, and phase distributions for the two simulated galactic conditions. New features of models enabled by the TIGRESS-NCR framework include maps of radiation fields and chemical abundances, as well as a full phase separation using both temperature and hydrogen abundances. Section 4 examines SFRs and the ISM state in the context of the pressureregulated, feedback-modulated (PRFM) star formation model (Ostriker & Kim 2022) by analyzing the midplane pressure components and their ratio to SFR surface density (feedback yields).
Section 5 summarizes our simulation results and discusses observational constraints, also situating our work within the context of recent star-forming ISM numerical studies.
In this section, we introduce the TIGRESS-NCR numerical framework. This is an extension of the original TIGRESS framework (Kim & Ostriker 2017, which we refer to as TIGRESS-classic hereafter) coupled with photochemistry and UV radiation transfer modules, as detailed in Kim et al. (2023). We begin by describing the governing equations (Section 2.1), and then briefly summarize the methods for treating star formation and SNe (Section 2.2), radiation transfer (Section 2.3), and photochemistry and cooling/heating (Section 2.4). Readers who are familiar with TIGRESS-classic can skip to the latter two sections for the new features.

Governing Equations
We solve the MHD equations in a local Cartesian rotating frame, with background galactic differential rotation treated via the so-called shearing-box approximation (Goldreich & Lynden-Bell 1965;Hawley et al. 1995). We use the grid-based code Athena to solve the ideal MHD equations Stone & Gardiner 2009), employing a high-order Godunov method combined with the constrained transport algorithm .
The conservation equations for gas mass, momentum, and total energy are, respectively, and ∂E tot ∂t The magnetic field evolution is governed by the induction equation without explicit resistivity (ideal MHD): while B must obey the divergence-free constraint In the above, ρ = µ H m H n H is the gas density, n H the number density of hydrogen nuclei, µ H the mean molecular weight per H nucleus, and m H the mass of a hydrogen atom; v and B are velocity and magnetic field vectors, respectively; P and P B = B · B/(8π) are thermal and magnetic pressure, respectively; E tot = ρv 2 /2+P/(γ −1)+P B is the total energy density, where γ = 5/3 is the adiabatic index. We explicitly follow non-equilibrium abundances of molecular (x H2 ) and ionized (x H + ) hydrogen by solving the transport of abundances with source terms, where ρ s = ρx s is the mass density of species s (H 2 or H + ), and C s is the net creation rate coefficient. On the right hand side of the momentum equation (Equation 2), we have source terms due to the total gravitational force (−ρ∇Φ), Coriolis force (−2ρv × Ω), and radiation force (f rad ) per unit volume. The total gravitational potential Φ = Φ sg + Φ ext (z) + Φ tidal (x) includes the self-gravitational potential obtained as the solution of Poisson's equation (including contributions from both gas and young star clusters, represented numerically as sink/star particles), the external gravitational potential in the vertical direction (fixed in time; see Kim & Ostriker (2017) for the exact form), and the tidal potential which gives rise to the differential rotation of the background flow (nonrigid body rotation); see below for the last. In the energy equation (Equation 3), we then have mechanical energy source terms associated with the gravity and radiation pressure forces (there is no work from Coriolis forces) in addition to the radiative heating and cooling terms (G − L). We solve Equation 7 using a Fast Fourier Transform method with shearing-periodic boundary conditions in the horizontal directions (Gammie 2001) and open boundary conditions in the vertical direction (Koyama & Ostriker 2009). We include newly formed stars' gravity using the particle mesh method by depositing the particle mass using a triangle-shaped cloud to obtain ρ sp (Gong & Ostriker 2013). The center of our domain corotates with the background galactic rotation speed at galactocentric radius R 0 , i.e., Ω = Ω 0ẑ , and we assume the galactic rotation curve is flat, i.e., the shear parameter q ≡ −d ln Ω/d ln R = 1. As a result, Φ tidal (x) = −qΩ 2 0 x 2 , where x is the local radial coordinate (x = 0 at the domain center). The source terms due to galactic differential rotation are included in the hyperbolic integrator by a semi-implicit method (Crank-Nicholson time differencing) as described by Stone & Gardiner (2010). The gravity source term is also included in the integrator, while the radiation force and cooling/heating source terms are included using an operator-split approach (see below).
The main new features in the TIGRESS-NCR framework are the explicit treatments of chemical processes and associated cooling and heating terms. This is in contrast to the TIGRESS-classic framework, in which the heating rate per hydrogen Γ is spatially constant (but time-variable), set via a simple scaling relative to the solar neighborhood value of the globally-attenuated instantaneous FUV radiation field as produced by star cluster particles (Kim et al. 2020a). In TIGRESSclassic, the cooling function Λ(T ) is only a function of temperature with a temperature dependent mean molecular weight µ(T ) combining Koyama & Inutsuka (2002) and Sutherland & Dopita (1993).
In this work, we directly calculate chemical reaction rates and cooling/heating rates from key microphysical processes that depend on gas properties (n H , T , v), species abundances (x s ), radiation fields in three UV bands (E rad ; see Section 2.3), and the CR ionization rate (ξ cr ). Explicitly, we have where Z g and Z d are the gas metallicity and dust abundance relative to solar neighborhood values. Details of these functions are provided in Kim et al. (2023). This paper assumes Z g = Z d = 1, corresponding to solar metallicity with abundances of Asplund et al. (2009) and fractional mass of metals Z g, = 0.014; and mass of grain material relative to gas 0.0081 (Weingartner & Draine 2001a). Although here we adopt globally constant values for Z g and Z d , in principle they can be determined locally (with appropriate treatments for metal enrichment and dust formation and destruction processes), and the TIGRESS-NCR framework is applicable for a wide range of Z g and Z d except for gas with very low metal and dust content in the early universe. 2 As detailed in Kim et al. (2023), we note that the heating and cooling functions that we adopt follow Wolfire et al. (1995Wolfire et al. ( , 2003 in all essential aspects and produce results consistent with theirs for solar neighborhood conditions, while also allowing for varying metal and dust abundances as well as UV and CR fluxes (see also Bialy & Sternberg 2019). Because our simulations are timedependent, they also make it possible to test the extent to which the predicted thermal equilibrium state is actually reached.

Star Formation and Supernovae
In addition to the source terms given on the right hand sides of Equations (1)-(3), we also include sink and source terms associated with star formation and SN feedback. The treatments of star formation and SN feedback using sink/star particles are identical to methods adopted for TIGRESS-classic.
We form and grow star cluster particles based on the sink particle treatment in Athena first introduced by Gong & Ostriker (2013) and modified further for the TIGRESS-classic framework (Kim et al. 2020a). Within the control volume (3 3 cells) surrounding a cell undergoing unresolved gravitational collapse, we create a sink particle by turning gas mass and momentum into a particle's mass and velocity. The collapse criteria are: (1) a cell's density is higher than a threshold Larson-Penston density depending on sound speed and numerical resolution, (2) the cell is at a local gravitational potential minimum, and (3) flows along all three Cartesian directions converge toward the cell. Each particle represents a star cluster (with typical mass > 10 3 M for our adopted resolution) consisting of coeval stars from a fully-sampled initial mass function (IMF). We use the STARBURST99 stellar population synthesis (SPS) model (Leitherer et al. 1999(Leitherer et al. , 2014 to determine SN rates for each star cluster, assuming a Kroupa (2001) IMF and Geneva evolutionary tracks for non-rotating stars.
Each star cluster hosts multiple SNe over its lifetime (t age < 40 Myr). We assume 50% of SN events are in binary OB systems, and if an event was from a binary we inject a massless particle with a velocity kick (Eldridge et al. 2011). These runaway stars can travel far from the cluster particle before they explode as SNe. However, we do not consider runaways as sources of UV radiation because the computational cost of ray tracing would become too expensive if runaway sources were included, and tests show that they do not contribute significantly to the total luminosity or ionization budget (Kado-Fong et al. 2020).
For each SN event, we first calculate the enclosed mass M SNR (sum of ejecta, M ej = 10 M , and preexisting ambient ISM) and mean density n amb of the surrounding ISM within a spherical region with a radius of 3∆x. If the calculated gas mass is less than the shell formation mass M sf = 1540 M (n amb / cm −3 ) −0.33 when a remnant becomes radiative (Kim & Ostriker 2015a), a total of E SN = 10 51 erg energy is injected with a thermal to kinetic energy ratio consistent with that of the Sedov-Taylor stage of evolution, after aver-aging out the feedback injection region. If the Sedov-Taylor stage is unresolved (i.e., M SNR > M sf ), a total of p SNR = 2.8 × 10 5 M km s −1 (n amb / cm −3 ) −0.17 radial momentum is injected. Given the high resolution and natural clustering of SNe realized in our simulations, only a small fraction of SN events (< 10%) are realized in the form of pure momentum injection.

UV Radiation Transfer and Cosmic Rays
All star cluster particles with (mass-weighted) age t age ≤ 20 Myr act as sources of UV radiation. Appendix C in Kim et al. (2023) provides details regarding radiation characteristics of the adopted SPS model.
We divide UV radiation into three frequency bins: photoelectric (PE; 110.8 nm < λ < 206.6 nm), Lyman-Werner (LW; 91.2 nm < λ < 110.8 nm), and Lyman continuum (LyC; λ < 91.2 nm). Both LW and PE photons (collectively referred to as FUV) provide an important source of gas heating via the photoelectric effect when absorbed by small dust grains and large molecules. All FUV photons are attenuated by dust along rays, and the LW band photons also dissociate H 2 and CO, and ionize C. To compute the dissociating radiation field for H 2 , we apply the Draine & Bertoldi (1996) self-shielding function to the LW band, using the H 2 column density integrated from each source to each cell. The LyC photons (also referred to as EUV) ionize neutral hydrogen (H and H 2 ) and are attenuated by dust.
To follow the propagation of UV photons from young star clusters, we utilize the adaptive ray tracing module in the Athena code (Kim et al. 2017b). After the hyperbolic part of the governing equations (including shearing-box and gravitational source terms) is integrated, we solve the time-independent UV radiation transfer equation,k for each frequency bin j along a set of rays. Here, I j is the intensity,k is the unit vector specifying the direction of radiation propagation, and α j is the absorption cross section per unit volume. In brief, 12×4 4 photon packets are injected for each band at the location of each source on a set of rays covering solid angles, corresponding to HEALPix level 4 (Górski et al. 2005). Photon packets propagate radially outward, and rays are split into four children when needed to ensure that each cell is intersected by at least 4 rays per source. The optical depth of each ray through each intersecting cell is computed and used to apply the corresponding rate of energy and momentum deposition. The radiation energy density E rad,j and flux F rad,j in each cell are obtained by summing over contributions from all rays passing through it. We then have f rad = j ρκ j F rad,j /c, which we use to update the gas momentum and corresponding kinetic energy density; the values of E rad,j in each cell are employed in photochemistry and heating calculations (Section 2.4). As with other fluid properties, the shearing-periodic boundary conditions are implemented for the ray tracing. Photon packets crossing the local-radial (x) edges of the computational domain are offset by the distance the boundaries have sheared in the local-azimuthal (ŷ) direction, and the position of sources is offset accordingly. The boundary condition in the y direction is periodic.
We terminate a ray if a photon packet exits the zboundary of the computational domain or the optical depth measured from the source is larger than τ max = 30 in all frequency bins. With just these basic ray termination conditions, however, we find that the computational cost of performing adaptive ray tracing every timestep is prohibitive. To reduce the cost of ray tracing without losing accuracy too much, we adopt three modifications.
(1) We perform ray tracing at intervals ∆t 2p based on the Courant-Friedrichs-Lewy (CFL) time step for the cold and warm gas at T < 2.0 × 10 4 K, or whenever a new radiation source is created, or whenever an existing radiation source is removed. (2) We put a hard limit on the maximum horizontal distance a ray may propagate from each source, which we denote as d xy,max . (3) We terminate a ray in the FUV band if the ratio between the luminosity of the photon packet and the total luminosity of all sources in the computational domain falls below a small number ε PP . The first condition, on the interval for radiation updates, is justified because the hot gas has very low opacity and its interaction with radiation is minimal. If there is no limitation on the maximum horizontal propagation distance of rays, the cost of ray tracing can become prohibitively high. We have found that imposing condition (2) reduces the cost to an acceptable level without significantly affecting the radiation field solution in the midplane region, provided d xy,max is large enough (see Section A.2). The condition (3) limits the maximum distance traveled by photon packets originating from faint sources, without seriously degrading the accuracy of the radiation field.
Terminating rays based on d xy,max and ε PP causes underestimation of the FUV radiation field at high altitudes. To address this issue without incurring additional computational cost, we apply an analytic solution based on the plane-parallel radiation transfer (see Section A.1). We stop ray-tracing for the PE and LW bands at |z| > z p-p and measure the area-averaged intensity of the PE and LW bands as a function of cos θ =k ·ẑ. We then calculate radiation energy density by integrating the intensity with the mean density averaged horizontally at each z. We adopt z p-p = 300 pc. This approach gives the mean radiation field as a function of z, which is uniform horizontally at a given z. It is generally adequate for high-z regions where the majority of gas is diffuse. For the LyC band, we do not apply the condition (3), nor do we adopt the plane-parallel approximation at large |z|.
The background CR ionization rate is scaled relative to the solar neighborhood level based on the SFR. Specifically, we adopt ξ cr,0 = 2×10 −16 s −1 Σ SFR,40 /Σ gas , where Σ SFR,40 and Σ gas are the SFR surface density measured from stars formed in last 40 Myr and instantaneous gas surface density relative to solar neighborhood values Σ SFR = 2.5 × 10 −3 M kpc −2 yr −1 and Σ gas = 10 M pc −2 . The local CR ionization rate is then set to where N eff is the effective shielding column density and N 0 = 9.35 × 10 20 cm −2 . To compute N eff at each zone, we additionally follow the radiation energy density in the PE band without dust attenuation. We then use the ratio of attenuated to unattenuated PE radiation energy density to obtain N eff = 10 21 cm −2 Z −1 d ln(E PE,unatt /E PE ). We note that CR transport in the ISM is uncertain and our prescription for CR attenuation in setting ξ cr should be considered provisional; the term 1/Σ gas in setting ξ cr,0 represents attenuation under average conditions and is motivated by Wolfire et al. (2003). The additional attenuation at high column in Equation (12) is motivated by observations of column-density dependent CR ionization rate (Neufeld & Wolfire 2017).

Photochemistry, Cooling, and Heating
In the MHD integrator, we transport molecular and ionized hydrogen using passive scalars (x H2 and x H + , respectively) without source terms as implemented in Athena. We obtain the atomic hydrogen abundance from the closure relation x H = 1 − 2x H2 − x H + . We then update the temperature and abundances in an operator split manner by solving two ordinary differential equations (ODEs): where e = P/(γ − 1) is the internal energy density, ) the mean mass per particle, and K = n H µ H k B /(γ − 1) is taken as a constant. Given the generally short cooling/heating and chemical time scales, in integrating these ODEs we take substeps (relative to the MHD time step) with the time step size determined by the minimum of MHD time step and 10% of the cooling, heating, and chemical time scales.
At each substep, we solve the two ODEs sequentially. Equation 13 is solved using the first-order backward difference formula with a Taylor expansion of the source terms, G − L, with respect to temperature which depends on the previous step's abundances and other information (see parameter dependence in Equation 9 and Equation 10 as well as Section 4 of Kim et al. (2023)). We then evaluate C s (see Section 3.1 of Kim et al. (2023)) and solve Equation 14 by treating it as a system of linear ODEs and use the backward Euler method. After the time-dependent update of hydrogen species, we compute the abundances of O + , C + , CO, C, and O under the steady-state assumption (see Sections 3.2 and 3.3 of Kim et al. (2023); see also Gong et al. 2017). We finally calculate the electron abundance x e from H + with contributions from C + , O + , and molecules at T < 2×10 4 K; or from He and metals assuming CIE at T > 2 × 10 4 K.
We refer the readers Kim et al. (2023) for complete information on the processes we include and rates we adopt. Here, we list the formation and destruction processes that are explicitly considered, as well as the cooling and heating processes.
• Ionized hydrogen: formation by photoionization, CR ionization, and collisional ionization of atomic hydrogen; destruction by radiative recombination and grain-assisted recombination.
• Ionized Carbon: formation by photoionization, CR-induced photionization, and CR ionization of atomic carbon; destruction by grain-assisted recombination, radiative+dielectronic recombination, and the CH + 2 formation reaction. • Heating: photoelectric effect on small grains by FUV photons; CR ionization of H and H 2 ; photoionization of H and H 2 ; formation, photodissociation, and UV pumping of H 2 .
H 2 ; bremsstrahlung and radiative/grain assisted recombination of free electrons with H + .
We smoothly transition from Λ others to Λ CIE at 2 × 10 4 K < T < 3.5 × 10 4 K using a sigmoid function, while Λ hyd is applied at all temperatures using time-dependent hydrogen abundances.
We note that for the dust-associated process, we adopt an empirically constrained dust population model of Weingartner & Draine (2001a), which consists of a separate population of carbonaceous and silicate grains as well as very small grains, including polycyclic aromatic hydrocarbon molecules. Our standard choice is their model with grain size distribution A, R V = 3.1, and b C = 4.0 × 10 −5 .

Models
We consider two galactic conditions R8 and LGR4, as described in Table 1, which are analogous to the models of the same names in the TIGRESS-classic suite (Kim et al. 2020a, the "LG" stands for the model with lower external gravity at a given gas surface density). The R8 model represents conditions similar to the solar neighborhood (e.g., McKee et al. 2015). In terms of gas and stellar surface densities, the conditions in models R8 and LGR4 roughly correspond to the areaweighted and molecular-mass-weighted averages of conditions in nearby star-forming galaxies surveyed as part of PHANGS (Sun et al. 2022).
For both simulations, the domain is a tall box, with dimensions (1024 pc) 2 × 6144 pc for R8, and (512 pc) 2 × 3072 pc for LGR4. We use d xy,max = 2048 pc and ε PP = 0 for R8 and d xy,max = 1024 pc and ε PP = 10 −8 for LGR4. With these choices of numerical parameters for raytracing, we found good convergence for the EUV radiation field everywhere, the FUV field near the midplane, and overall simulation outcomes (see Section A.2). We note that the selection of the optimal values of d xy,max and ε PP depends on the system's conditions, especially on Z d (which determines dust absorption).
We run each simulation with two different resolutions: 8 pc and 4 pc for R8; and 4 pc and 2 pc for LGR4. The initial gas distribution follows double Gaussian profiles (see Kim & Ostriker 2017) representing warm and hot components with the Gaussian scale height corresponding to initial velocity dispersions of 10 and 100 km/s for R8 and 15 and 150 km/s for LGR4. To reduce initial transients, we add additional velocity perturbations with amplitude of 10 and 15 km/s for R8 and LGR4, respectively, along with randomly distributed initial star clusters that provide non-zero initial heating and SNe. The initial magnetic field is set to be azimuthal B = B 0ŷ with uniform plasma beta β 0 ≡ P th /(B 2 0 /8π) = 1, which is close to the expected saturation value of the magnetic field (e.g. Kim & Ostriker 2015b;Ostriker & Kim 2022). After one or two cycles of star formation and feedback, the system reaches a quasi-steady state with minimal impact from initial transients and the choice of density profiles and the level of initial turbulence.

TIGRESS-NCR SIMULATIONS
We first provide an overview of the R8 and LGR4 simulations. Here, we focus on global ISM properties and energetics as well as the distribution of the multiphase (both thermally and chemically) ISM near the galactic midplane. Figure 1 shows the time evolution of key global quantities including (a) the gas surface density Σ gas ≡ M gas /(L x L y ) along with surface densities of newly formed stars Σ * ,new ≡ M * ,new /(L x L y ) and mass loss via outflows from the computational domain Σ out ≡ M out /(L x L y ); (b) the SFR surface density over the last ∆t = 10 Myr

Global ISM Properties
LGR4-2pc 37.9 +1. where M * ,i (t age < ∆t) is the total mass of star particles with age younger than ∆t; (c) the effective vertical velocity dispersion Here, the values of σ eff and H are calculated only for the warm and cold gas with temperature T < 3.5 × 10 4 K. We discuss phase separated velocity dispersions in Section 3.4. We note that the magnetic term in σ eff is not simply the magnetic pressure, but instead represents the vertical component of Maxwell stress including both magnetic pressure and tension. We measure the mass-weighted vertical velocity dispersion of the turbulence for the warm and cold gas, i.e. σ z,turb ≡ ( ρv 2 z dV / ρdV ) 1/2 , and report this separately in Table 2.
The simulations begin with initial rough hydrostatic equilibrium with H ∼ 150 − 200 pc. The idealized initial setups soon lead to a burst of star formation and associated feedback as the disk loses its initial vertical support from both thermal and turbulent pressure. During the first ∼ 100 Myr evolution, both models experience at least two complete star formation and feedback cycles (rise and fall of SFRs), whose timescale is proportional to the vertical crossing timescale of the disk (Kim et al. 2020a). To reduce the computational time needed for high-resolution models, we refine and restart low-resolution simulations (R8-8pc and LGR4-4pc) after the initial transient has passed (t = 200 Myr). The meshrefined, restarted models are run for longer than one orbit time (or 3-4 star formation and feedback cycles) to obtain a fair statistical sample of states at higher resolution. In Table 2, we summarize mean and standard deviation of quantities of interest over t = 250−450 Myr and t = 250 − 350 Myr for R8 and LGR4, respectively, at different resolutions. Our results verify the overall convergence of the global properties with respect to numerical resolution. 3 As shown in the top row of Figure 1, the gas surface density (solid) decreases gradually due to star formation (dashed) and outflows (dotted). The global properties shown in the bottom three rows of Figure 1 reach a quasi-steady state, with substantial temporal fluctuations (∼ 0.2 − 0.3 dex), and show quasi-periodic fluctuations. The characteristic period is the vertical oscillation time determined by the total (gas+star+dark matter) midplane density ∼ (4πGρ tot ) −1/2 , which is similar to the vertical crossing time (see Kim et al. 2020a). In Table 2, we list the vertical crossing time t ver ≡ H/σ z,turb and gas depletion time t dep ≡ Σ gas /Σ SFR . Quasiperiodic fluctuating behavior in Σ SFR and σ eff shows higher frequency fluctuations than H. Occasionally, when there is a big burst, systematic offsets among three quantities stand out; a peak of SFR is followed by a peak of velocity dispersion and then scale height (e.g., see peaks near 100 Myr and 300 Myr in R8-8pc).

Global Energetics
Figure 2(a) shows the total energy injection rate per unit area by UV radiation 4Ṡ UV ≡ L UV,tot /(L x L y ) as a function of time, where L UV,tot is the total UV (PE+LW+LyC) luminosity of star particles with t age < 20 Myr. This energy injection rate is determined by the adopted SPS model (STARBURST99 in our case) and recent star formation history. About 5% of UV radiation energy goes into heating the warm and cold gas. The total radiative cooling always exceeds radiative heating because the cooling offsets heating provided by SNe. Only 2-3% of the injected SN energy is advected out of our simulation domain.
UV radiation propagates through the ISM and is absorbed by gas and dust, photoionizing some regions where EUV penetrates, and heating up the gas via the photoelectric effect in other regions where FUV penetrates. In addition, CR ionization is an important heating source in regions that are shielded to both EUV and FUV. The total radiative (including CR) heating rate per unit areaṠ heat ≡ GdV /(L x L y ) shown in Figure 2(b) is the sum of hydrogen photoionization heating by LyC radiation (Ṡ heat,PI /Ṡ heat ∼ 75%), photoelectric heating from FUV (PE+LW) radiation on small dust grains (Ṡ heat,PE /Ṡ heat ∼ 20%), and CR ionization heating (Ṡ heat,CR /Ṡ heat ∼ 1 − 2%), plus a tiny contribution from H 2 heating (< 0.1%). The global heating efficiency, defined as the ratio of the total heating rate to the UV radiation injection rate, isṠ heat, The radiative heating within the simulation domain is balanced by radiative cooling. Figure 2(c) shows the difference between cooling and heating per unit area within the simulation volume. The difference is positive, indicating net cooling. The excess in radiative cooling is offset by energy input from SN feedback (Ṡ SN ; Figure 2(d)).Ṡ SN is about two orders of magnitude smaller than the UV radiation injected (Ṡ UV ; Figure 2(a)), and a factor ∼ 3 lower than the radiative heating rateṠ heat . A small fraction of energy also leaves the computational domain through outflows; the majority of outflowing energy is in the hot gas and therefore originally deposited by SNe, and the kinetic energy of outflowing gas is also powered by SNe. This outflowing energy escaping the domain (∼ 2 − 3% ofṠ SN ) accounts for the small excess of SN injection energy over the net cooling within the domain.

ISM Cartography
The instantaneous spatial distribution of the ISM mass and energy densities is highly structured and complex. To provide a visual impression of the ISM structure in our simulations, we display maps of various quantities from R8-4pc at two epochs, t = 280 and 295 Myr in Figure 3(a) and (b). These times respectively correspond to a local peak and trough of Σ SFR (see Figure 1(c)). We note that qualitative features presented here for R8 are also seen in LGR4.
We first focus on the epoch shown in Figure 3(a), shortly after a burst of star formation has occurred, during which many new star clusters were formed. Very young clusters (t age < 5 Myr) act as strong UV radiation sources. These clusters are spatially correlated with each other and with the dense clouds where they were born. There are two distinct types of bubble structures: hot SN-driven bubbles (characterized by low density, diffuse EM, and high temperature) and warm ionized bubbles (characterized by high EM). The electron fraction of the two types of bubbles are different. The hot ionized gas has higher x e ≈ 1.2 (bright green), due to free electrons from collisionally ionized H, He, and metals. In contrast, x e ≈ 1 (dark green) in the warm ionized gas as electrons are mostly from photoionized H (and   (bottom) in the x-y plane (corresponding to integrals of ρ and n 2 e along the z-axis, respectively). The letters in the emission measure map indicate regions of ionized bubbles (warm ionized bubble -A1, A2, A3; young superbubble -B1, B2; old superbubble -C). Right: Slices through the midplane, z = 0. From left to right, the top row shows hydrogen number density nH, temperature T , and electron fraction xe; the bottom row shows normalized FUV radiation energy density χFUV, cooling rate L, and net heating rate G − L. χFUV ≡ JFUV/J D78 FUV , where JFUV is the FUV mean intensity and J D78 FUV = 2.1 × 10 −4 erg s −1 cm −2 sr −1 (Draine 1978). Note that for G − L, the pink colormap is used only for positive values (net heating), while the blue colormap from L is used for negative values (net cooling). Contours of EUV radiation energy density are overlaid in the χFUV panels for log10[ELyC/( erg cm −3 )] = −15 (red), −14 (orange), −13 (green), and −12 (blue). In the Σgas panel, young star clusters with age < 40 Myr and |z| < 50 pc are shown as circles. The size of circles is proportional to the cluster mass, but in practice its range is narrow (10 3 − 5 × 10 3 M ). Clusters with age < 20 Myr (magenta-ish colors; see colormap in the bottom-right of the Σgas panel) are FUV sources, while very young clusters (age < 5 Myr) are the only significant EUV sources (enclosed by green/blue contours in the χFUV panel). Locations where SNe exploded within the past 10 Myr, and within |z| < 15 pc of the slice shown, are marked as star symbols in the nH panel.
a tiny contribution from C + , O + , and molecular ions).
In the upper region of the domain (y ∼ 0.2 kpc), examples of warm ionized bubbles, corresponding to high-EM sites, are marked as A1, A2, and A3. In the middle region (y ∼ 0 kpc), two superbubbles that are formed relatively recently and show moderate EM are marked as B1 and B2. An example of an old, low EM superbubble is marked as C.
It is evident that recently born star cluster complexes are responsible for photoionizing bubbles A1, A2, and A3. Ionizing radiation from clusters near A1 and A2 is fairly well blocked toward bubbles B1 and B2, while extended radiation from these sources could ionize a substantial area toward the top of the domain. Large areas within the domain remain neutral (x e 10%; whiteblue in the x e panel) as EUV is effectively truncated due to the large cross section of the neutral hydrogen. FUV is only significantly absorbed by dense clouds, making the radiation energy density low in their interiors and also casting shadows behind them. Still, it is evident that most of the neutral gas is irradiated with FUV (χ FUV 0.5), which is a major heating source. Bubble C is an old, hot superbubble, and star clusters are old and no longer contributing significant EUV. As a result, the hot gas is bounded by the neutral gas. In contrast, the intermediate age clusters inside hot bubbles B1 and B2, together with other nearby young clusters, create a photoionized region between the hot and neutral gas. The strongest cooling in the slice, as shown in the L panel, occurs in the photoionized regions near Bubbles A and B; the main coolants are nebular lines of metal ions. However, the heating produced by ionizing radiation offsets or even exceeds the cooling in this region, leading to net heating (see the G − L panel). The net cooling rate is highest at hot bubble boundaries (CIE cooling at T ∼ 10 5 K). These interfaces where hot gas mixes with denser gas to become strongly radiative are important in bubble energetics.
As young star clusters shown in Figure 3(a) age, they begin to produce SNe, resulting in superbubbles, which merge into a very large hot bubble in the center of Figure 3(b). This is carved by several clustered SNe (positions are shown as star symbols in the n H panel). At this epoch, there is only one significant ionizing source at the center of the midplane slice, and the area filled with the warm ionized gas (dark green in the x e panel) is greatly reduced. There are a few out-ofmidplane sources (not shown in the Σ gas panel) responsible for large EM bubbles at (x, y) ∼ (0, 0.5) kpc and (x, y) ∼ (−0.2, −0.3) kpc. It is also noticeable that old clusters are now spread across the simulation domain; clustering of clusters is reduced.
There is temporal evolution in the ISM phase structure over the interval shown between the two snapshots. The midplane volume filling factor of the warm ionized medium achieves its local maximum, ∼ 30%, at t = 280 Myr (Figure 3(a)), decreasing to ∼ 10% within another 5 Myr. In contrast, the midplane hot gas filling factor increases gradually from 20% at t = 280 Myr to 50% at t = 295 Myr. The filling factor of the warm neutral medium changes from 40% to 20% over the same 15 Myr interval.

Phase Definition, Filling Factor, and Velocity Dispersion
Traditionally, in the ISM literature, the gas is often divided into different phases based on temperature and hydrogen chemical state. We choose a set of specific temperature and abundance cuts to define 9 phases as summarized in Table 3. The distributions of these phases are depicted for a sample snapshot from LGR4-2pc in Figure 4. In our previous work (Kim & Ostriker 2017 and subsequent works based on TIGRESS-classic), we used temperature cuts only to define 5 ISM phases. Key additional information available in the current study is the time-dependent hydrogen abundance, allowing for subdivisions, e.g., "warm" into the warm neutral and warm ionized medium. The warm ionized medium is further divided into "warm-photoionized" and "warmcollisionally-ionized" medium with a temperature cut at T = 1.5 · 10 4 K, above which hydrogen can be collisionally ionized (see red dashed line in Figure 4(c)). We assign every cell to one of the 9 phases exclusively.
A summary of the mass and volume fraction for each phase is shown in Table 4, for both R8 and LGR4. Here, we do not explicitly distinguish CNM and CMM and we ignore UIM given its negligible mass and volume fractions. Instead, we use Cold for the combined cold medium (CMM+CNM). We note that, as shown in Figure 4(c), hydrogen species abundances vary continuously, and a significant amount of partially molecular gas is present in CNM. The total molecular gas mass is thus larger than the mass of CMM. We note that the Cold mass fraction increases substantially with higher numerical resolution at the expense of UNM and WNM, but the fractions of all the other phases are reasonably converged. WNM fills the majority of the volume, with substantial contributions from WPIM and HIM as well as UNM. The neutral gas (Cold, UNM, and WNM) dominates the total mass budget. WPIM contributes to both volume and mass at ∼ 10% level, with an increasing contribution at high altitudes (e.g., Kado-Fong et al. 2020;see also N. Linzer et al. in prep).
b In principle, 'neutral' includes both 'atomic' and 'molecular'. But historically, the cold neutral medium has been used to denote cold atomic medium. Here, we follow the convention. c This includes cold temperature range but is dominated by unstable. Regions assigned to mutually exclusive defined phases as shown in the key. (c) Mass-weighted joint PDF of log10 T and xH for gas within |z| < 300 pc, with dividing lines for different gas phase definitions (see Table 3). The red dashed curve in panel (c) denotes the relation between xH and T based on CIE of hydrogen. effective vertical velocity dispersion as defined by Equation 16 and the mass-weighted turbulent velocity dispersion only considering the P turb = ρv 2 z term for each component. Given that the neutral medium (especially, warmer component UNM+WNM) dominates the mass fraction (Table 4), σ eff,2p agrees well with the effective velocity dispersion of all warm and cold gas at T < 3.5×10 4 K presented in Table 2. On the one hand, this makes the observed H I velocity dispersion a good tracer for the (thermal plus turbulent) velocity dispersion of the massdominating component. On the other hand, it shows that WIM tracers will typically overestimate the massweighted velocity dispersion. This could lead to a bias if, for example, Hα velocities are used in estimators for the ISM weight (see Section 4). It is also noteworthy that the turbulent velocity dispersion is much lower than the effective velocity dispersion; thermal and magnetic components contribute significantly to the total pressure. Figure 5 and Figure 6 show probability distribution functions (PDFs) of temperature, density, and thermal pressure from R8-4pc and LGR4-2pc, respectively, based on the region within |z| < 300 pc. WNM and WPIM have LGR4-4pc 2.1 +1.0  similar characteristic densities and temperatures, but the thermal pressure of WPIM is higher because of the contribution from free electrons. CMM corresponds to the dense part of CNM with similar thermal pressure. All other phase definitions are mostly equivalent to simple temperature cuts. UNM and WNM are in rough thermal pressure equilibrium, but CNM tends to have lower thermal pressure, which is compensated by higher magnetic pressure. WCIM and WHIM are not significant components in terms of mass and volume as they are usually populated only near the interfaces between warm and hot gas (see Figure 4). However, most (net) cooling occurs in these phases (see bottom right panel of Figure 3, and Section 3.5 below). The thermal pressure of HIM is generally larger than that of WNM. Since the thermal pressure of HIM is in balance with the total pressure of WNM, and the turbulent and magnetic contributions in WNM are larger in LGR4 than in R8 (see Section 4), the difference in thermal pressure between HIM and WNM is larger in LGR4.
3.5. Joint PDF of density and pressure Figure 7 shows, for R8-4pc at t = 320 Myr within |z| < 300 pc, the instantaneous distribution of gas in the density-pressure phase plane weighted by volume, mass, and net cooling rate. The total gas is shown in column (a), while column (b) shows the neutral (atomic + molecular) gas and (c) shows the ionized gas, using a cut x H + = 0.5. Note that the x-axis is the hydrogen number density n H rather than the total number density n = (1.1 + x e − 0.5x H2 )n H . Therefore, at a given temperature the neutral and ionized medium lie on different pressure tracks as a function of n H -a higher (lower) pressure track for the warm ionized (neutral) gas.
In TIGRESS-NCR, the heating and cooling rates are not solely a function of density and temperature and a spatially-uniform FUV field (which was the case in TIGRESS-classic), but also depend on other quantities such as the electron abundance and spatiallynonuniform radiation (see Equation 9 and Equation 10).    Figure 3(a)). All gas within |z| < 300 pc shown in (a) is subdivided into (b) neutral (x H + < 0.5) and (c) ionized (x H + > 0.5). Note that the PDF weighted by net cooling rate is normalized by the total cooling rate within the volume, adopting logarithmic blue-ish and pink-ish color scales for net cooling and heating, respectively. In the middle column, the diagonal dotted lines show T = 500 and 6000 K for neutral gas (P/kB = 1.1nHT ). In the right column, the diagonal dotted lines show T = 3.5 · 10 4 and T = 5 · 10 5 K for ionized gas (P/kB = (1.1 + xe)nHT with xe = 1 and 1.2, respectively). The majority of the WPIM is found near T ∼ 7000 K. The neutral medium is distributed somewhat broadly, but with a concentration near T = 7500 K for WNM, and near the locus corresponding to thermal equilibrium (zero net cooling) for CNM. The HIM region shows tracks roughly following P ∝ ρ 5/3 , corresponding to adiabatic expansion of hot bubble interiors. Horizontal dashed lines in (b) and (c) show volume-weighted mean pressure of the neutral gas and ionized gas, respectively. In the middle column (b), we show unshielded equilibrium curves for ξcr = 2.9 × 10 −16 s −1 and three values of FUV radiation field χFUV = 0.51, 0.87, and 2.6, corresponding to 2 nd , 50 th , and 98 th percentiles of the volume-weighted χFUV distribution within |z| < 300 pc. The median curve describes the WNM branch well.
Thus, a single thermal equilibrium curve applicable for the whole simulation domain cannot be drawn in Figure 7. Yet, we still see the characteristic locus of thermal equilibrium for neutral gas (see e.g., Field et al. 1969;Wolfire et al. 1995;Kim et al. 2023) in the bottom panel of Figure 7(b) as the boundary between coolingdominated and heating-dominated regions. Given ξ cr = 2.9 × 10 −16 s −1 for the background gas and the median value of χ FUV = 0.87, in Figure 7(b), we plot an equilibrium curve as thick solid line (as well as two thin lines for χ FUV = 0.51 and 2.6). Since we ignore shielding of FUV in these one zone models, the unshielded equilibrium curves give overall higher pressure at high densities, although the WNM equilibrium branch and its maximum pressure is well described by the median equilibrium curve.
The volume-weighted mean pressure (within |z| < 300 pc) of the neutral gas, P/k B = 3.1 × 10 3 cm −3 K, is shown as a horizontal dashed line. This pressure sits nicely between the maximum thermal equilibrium pressure of WNM (the bulk net heating region shown as pink) at n H ∼ 0.5 cm −3 and P/k B ∼ 5 × 10 3 cm −3 K and the minimum thermal equilibrium pressure of the CNM (the bulk net cooling region shown as blue) at n H ∼ 5 cm −3 and P/k B ∼ 1 × 10 3 cm −3 K. Although the ionized gas has a very wide range of thermal pressure, the mean is P/k B = 7.9 × 10 3 cm −3 K, which is shown as the horizontal dashed line in Figure 7(c). This is higher than that of the neutral gas, in which turbulence and magnetic field significantly contribute to the total pressure (see Section 4.1).
As shown in Figure 5, both mass and volume are dominated by the neutral medium near the disk midplane. The ionized gas occupies ∼ 30 − 40% by volume (approximately equally in WIM and HIM) and ∼ 10% by mass (mostly in WIM). The bottom row of Figure 7 shows that both neutral and ionized gas populate a wide parameter space with net cooling or heating (i.e, gas is out of thermal equilibrium). Photoionization can cause net heating in expanding H II regions as well as the diffuse WIM, which is evident in the narrow dark pink strip at T ∼ 7000 K in the bottom-right panel of Figure 7. Out-of-equilibrium CNM is mostly in the net (radiative) cooling regime, implying that dissipation of turbulence may contribute to the thermal balance in CNM to allow an overall excess of radiative cooling over radiative heating. Locally, WNM is perturbed into both net cooling and heating regimes, while WNM as a whole is experiencing net cooling.
Due to its low density, cooling in the hot gas located inside SN(e)-driven bubbles is negligible. The corresponding high-temperature regions in Figure 7 show adi-abatic expansion tracks, P ∝ ρ 5/3 . The thermalized SN energy is mostly cooled away at lower temperature phases, e.g., WHIM and WNM/WIM. To clarify the contribution of each gas phase in cooling, Figure 8 plots the cooling and net cooling contribution within |z| < 300 pc from each phase as a function of density, for R8-4pc (left) and LGR4-2pc (right). The total radiative cooling rate per volume, shown in the top panel, is dominated by WPIM (∼ 70%). WNM and WHIM contribute about 10% each. But, WPIM and WNM are also the gas phases within which most radiative (photoionization and photoelectric) heating occurs. The bottom panels show the net cooling rate per volume, with heating subtracted from cooling. The net cooling, which when integrated over volume produces the history shown in Figure 2(c), is now dominated by WHIM (∼ 40%). WNM, WPIM, and WCIM contribute about 10-20% each.

PRESSURE-REGULATED, FEEDBACK-MODULATED THEORY OF THE EQUILIBRIUM STAR-FORMING ISM
Having presented the overall characteristics of the simulated ISM, in this section we focus on the midplane pressure and stresses, gas weight, and SFR surface density, and their mutual correlations. These analyses aim to confirm the validity and prediction of the PRFM star formation theory, first formulated in Ostriker et al. (2010) and Ostriker & Shetty (2011) and tested in subsequent work (Kim et al. 2011(Kim et al. , 2013Shetty & Ostriker 2012;Kim & Ostriker 2015b).
For a comprenhensive summary and detailed derivation of the theory, the reader is referred to Ostriker & Kim (2022). We closely follow the analysis of Ostriker & Kim (2022), which analyzes the TIGRESS-classic simulation suite.
The PRFM theory views the ISM in galactic disks as a long-lived thermal-dynamical system with stellar feedback as the main energy source. Despite the vigorous dynamical evolution seen in our simulations (and in the real ISM), the system is in a quasi-steady state on average (in terms either of long-term temporal averages or large-scale ensemble averages). Under this assumption, the governing gas dynamics equations dictate vertical dynamical equilibrium, a balance between total pressure and gas weight. The ISM energy would drop quickly through cooling and dissipation in the absence of inputs, but stellar feedback can offset losses to maintain the required pressure/stress. As a consequence, the PRFM theory posits that galactic SFRs are naturally linked to the dynamical equilibrium pressure, which in turn predicts galactic SFRs from large-scale mean galactic parameters. Numerical simulations that directly capture selfconsistent energy injection (by feedback) and loss processes are critical in determining the feedback yield. In the TIGRESS-NCR framework, we do not impose radiation fields and the resulting photoheating based on observational estimates, but compute J FUV via ray-tracing from young star cluster particles formed in our simulations, where the number and location of these star particles is self-consistently set by the rate of gravitational collapse. Similarly, our simulations have sufficient resolution to follow the transition from adiabatic to cooling stages of SN remnant expansion. Thus, unlike in lower resolution simulations, we resolve the simultaneous heating and acceleration of ambient gas by SN shocks as well as subsequent cooling. The present simulations solving direct UV radiation transfer not only for FUV but also for EUV are critical for validation of the simpler treat-ment for FUV heating used in TIGRESS-classic, and for accurately calibrating the feedback yields.
Our analysis steps in this section are as follows. We first check pressure equilibrium among the three phases and the vertical dynamical equilibrium between total pressure support and gas weight (Section 4.1). Then, we measure each pressure component and examine the pressure-Σ SFR relation (Section 4.2). This gives a numerical calibration of the key parameter in the theory, the ratio of pressure and Σ SFR , which we call the feedback yield. We compare our new results for the feedback yield with the theoretical and numerical results in Ostriker & Kim (2022). Since we only consider two galactic conditions in this paper, we refrain from deriving new fitting results. In Section 4.3, we present the correlations between SFR surface density, total pressure, and weight (or its simplified estimator, dynamical equilibrium pressure).

Pressure Equilibrium and Vertical Dynamical Equilibrium
By integrating the vertical component of the momentum equation (Equation 2) from the midplane to the top/bottom of the gas disk, the vertical dynamical equilibrium condition assuming a steady state is then given by the balance between the midplane total pressure (P tot ) and the ISM weight (W): ∆ P tot ≡ ∆ P th + P turb + Π mag + ∆P rad = W (18) where ∆ denotes the difference between the midplane (z = 0) and the top/bottom of the gaseous disk (z = ±L z /2), and the angle brackets denote a horizontal average. In Equation 18, we adopt the following nomenclature of pressure components: thermal pressure P th = P , turbulent pressure (Reynolds stress) P turb ≡ ρv 2 z , and magnetic stress (magnetic pressure + tension) Π mag ≡ (B 2 x + B 2 y − B 2 z )/(8π). Note that the mean or turbulent magnetic stress (Π B or Π δB ) is respectively defined by substituting for B with the mean B ≡ B or turbulent δB ≡ B − B component. The radiation pressure and weight terms can be defined toward either upper or lower disk, ∆P rad = ±Lz/2 0 f rad,z dz and W = ±Lz/2 0 ρ∂Φ/∂z dz. We take an average of two values (integrated from top or bottom) in the following analysis. The vertical gravity −∂Φ/∂z is a sum of terms from stars, dark matter, and gas, so the total weight can be decomposed into two terms: gas weight from external gravity W ext (due to stellar disk and dark matter halo), and gas weight from self-gravity W sg (due to gas). Generally, the pressure components at the midplane z = 0 are much larger than those at the top of the gaseous disk, leading to ∆P → P (z = 0).
To separate the contribution from each phase, we define the horizontal average of a quantity q for a selected phase by where Θ(ph) = 1 if temperature and abundance condition in Table 3 is satisfied and 0 otherwise. Here, we simplify our full phase definition into three phases: 2p for the neutral medium at cold to warm temperatures (i.e., 2p=CMM+CNM+UNM+WNM), WIM for the warm ionized medium (i.e., WIM=WPIM+WCIM), and Hot for the hot ionized medium (i.e., Hot=WHIM+HIM). We can also separately define the area filling factor f A,ph ≡ Θ(ph)dxdy/L x L y . Each phase's contribution adds up such that P tot = ph P tot ph . Individual pressure components ( P th , etc.) can be written as a sum over contributions from each phase in the same way. The typical midplane value of the total pressure in each phase is defined bỹ P tot,ph ≡ P tot ph /f A,ph (and similarly forP th,ph etc.). We can then write P tot = ph f A,phPtot,ph . If the typical values of the total pressure for 2p, WIM, and Hot are comparable with each other, we have P tot ≈ P tot,X ph f A,ph =P tot,X , where X denotes any given phase. If we neglect the direct UV radiation pressure ∆P rad for succinctness as it contributes less than 1% to the total pressure in both models, Equation 18 simply becomes ∆ P tot ≈P tot,2p =P th,2p +P turb,2p +Π mag,2p = W.
(20) We note that weight (and radiation pressure) is vertically integrated over all phases, and that the (timeaveraged) pressure of gas in any given phase at the midplane must support the weight of gas in all phases above it (rather than selectively supporting its own phase). In our simulations, the gas weight is dominated by 2p with 14% (4%) from WIM for R8 (LGR4) and less than 1% from Hot. The contribution from the external gravity is 75% for R8 and 30% for LGR4. Figure 9 shows time evolution of all pressure terms in Equation 20 along with the total midplane pressures of WIM and Hot. In this and other figures and tables, the values of pressures shown represent midplane averages either within a given phase or over all phases, dropping the tilde for cleaner notation. Comparing the total pressures of each phase (dark blue for 2p, yellow for WIM, and red for Hot), we confirm that they are roughly in pressure equilibrium. Also shown are the direct measure of the ISM weight (W) and the commonly used weight estimator (which assumes that the stellar disk is thicker than the gaseous disk; Ostriker & Kim 2022) constructed from observables (e.g., Sun et al. 2020). We note that σ eff,2p presented in Table 5 is the massweighted mean for the 2p phase over the entire domain (not the midplane measure). This kinetic (ther-mal+turbulent) velocity dispersion is a direct observable given line emission from the neutral (atomic and molecular) gas, and then σ eff,2p can be obtained by correcting for the magnetic terms. On average, the total pressure and weight are in good agreement with each other (they are usually off-phased). Figure 10 plots P tot,2p as a function of (a) P tot,hot , (b) W, and (c) P DE,2p , while Table 6 summarizes the midplane pressure components in 2p as well as total pressure in each phase, weight, and weight estimator. Again, approximate pressure equilibrium among the different phases holds, but Hot gives slightly lower total The simple weight estimator PDE (grey dashed) provides reasonable agreement with the weight and total pressure. We show each pressure component of 2p as blue (P th,2p ), orange (P turb,2p ), and green (Πmag,2p) lines. LGR4-2pc 2.1 +0.4 pressure. Thermal pressure dominates in WIM and Hot, while thermal pressure is the smallest component in 2p, with magnetic and turbulent components comparable. P tot,2p ≈ W directly demonstrates that the ISM pressure is regulated in disk systems as it obeys the conservation "law" of momentum (on average). Figure 10(c) demonstrates the validity of P DE,2p as a reasonable estimator of the true weight (see Table 6) and hence total midplane pressure.

Feedback Modulation and Yields
The PRFM theory postulates that thermal and turbulent pressure (∝ energy density) components are sourced by feedback from massive young stars through heating by UV radiation and turbulence driven by SNe. The balance between radiative cooling and heating sets the thermal pressure, while the balance between turbulence driv-ing and dissipation sets the turbulent pressure. Magnetic fields are set by the saturation of galactic dynamo, providing the contribution from the magnetic term at a level similar to (or slightly below) the turbulent term (Kim & Ostriker 2015b). The pressure components are thus expected to scale with the rate of feedback energy injection, and therefore with Σ SFR . We define the feedback yields Υ c ≡ P c /Σ SFR as the ratios of a pressure component c to Σ SFR , to quantify the feedback modulation of individual pressure components. Note that the natural unit for the feedback yield is velocity.

Thermal Pressure
Because of the short cooling time in the cold and warm ISM (2p and WIM), the energy gains from radiative heating are quickly lost through optically thin radiative cooling mostly within the same phase. Thermal pressure in both 2p and WIM is then expected to scale with the radiative heating rate, which is proportional to the mean UV intensity and hence SFR. In the 2p medium, the photoelectric effect by FUV is the main heating source. Therefore, P th,2p ∝ Γ PE ∝ PE J FUV , where PE is the photoelectric heating efficiency, which depends sensitively on the grain size distribution and the grain charging parameter ψ ≡ G 0 √ T /n e (e.g., Bakes & Tielens 1994;Weingartner & Draine 2001b). The source of FUV radiation is massive young stars (the luminosityweighted mean age of FUV emitters is ∼ 10 Myr) so that J FUV = f τṠFUV /(4π) for f τ a factor accounting for UV radiative transfer in the dusty ISM, wherė A simple radiation transfer solution for uniformlydistributed sources in a uniform slab gives (Ostriker et al. 2010). Here, E 2 is the second exponential integral and τ ⊥ = κ FUV Σ gas is the mean optical depth to FUV. For κ FUV = 10 3 cm 2 g −1 , f τ ≈ 1/τ ⊥ at Σ gas > 20 M pc −2 . In the TIGRESS-classic suite, we adopted the approximate form of f τ as presented in Equation 22 to convertṠ FUV to J FUV , and we also adopted a single value for PE . The attenuation of FUV increases at higher surface densities (which generally corresponds to higher pressures). The relationship between the thermal pressure 5 We note that we have changed notation for the FUV luminosity per area from Σ FUV (e.g., Ostriker et al. 2010;Ostriker & Kim 2022) toṠ FUV to consistently refer to areal energy gain and loss rates usingṠ (see Section 3.2) and Σ SFR is thus sublinear, resulting in a decrease of thermal pressure yield at higher P DE . The fit to the TIGRESS-classic suite gives (Ostriker & Kim 2022) log 10 (P th,2p /k B ) = 0.603 log 10 (Σ SFR ) + 4.99 (23a) log 10 Υ th = −0.506 log 10 (P DE /k B ) + 4.45. (23b) In the current simulations, the connection from Σ SFR to J FUV to Γ PE is self-consistently determined by explicit UV radiation transfer and an adopted theoretical dust model for the heating efficiency and cross sections. Figure 11(a) plots P th,2p vs. Σ SFR , showing the similar sublinear relationship calibrated to TIGRESSclassic (Equation 23a; solid black line). Figure 11(d) shows the thermal feedback yield that decreases as P DE increases (the TIGRESS-classic fit Equation 23b is also shown). We find that the scaling is quite similar to the fit from the TIGRESS-classic suite, but the normalization in Υ th is higher here -that is, the TIGRESS-NCR simulations give rise to the higher thermal pressure at a given Σ SFR than the TIGRESS-classic simulations. The offset is because the explicit treatment of the heating efficiency here yields on average a factor of 2-3 higher heating rate for a given J FUV , compared to the heating rate coefficient adopted in TIGRESS-classic (which is from Koyama & Inutsuka 2002). The consistent scaling suggests that Equation 22 is indeed a good approximation for the mean attenuation factor in comparison to the actual radiation transfer solution obtained here by adaptive ray tracing (N. Linzer et al. in prep.).

Turbulent Pressure
The turbulent pressure in the warm and cold components of the ISM arises from large-scale forcing, with expanding hot bubbles produced by SN feedback as the most important source. Because the energy injection from SNe is highly localized in space and time, it creates a shock when it is transferred to the warm and cold ISM gas. This accelerates the surrounding ISM, increasing the total momentum until the shock becomes radiative when the post-shock temperature T 10 6 K (or v SNR ∼ 200 km s −1 ). The radial momentum injected per SN (p * ) is much larger (∼ 10 5 M km s −1 ) than the momentum of the initial SN ejecta (∼ 10 4 M km s −1 ) because the shock accelerates two orders of magnitude more mass than the initial ejecta before becoming radiative. The SN momentum injection is also much greater than other sources, such as expanding H II regions ) and stellar wind driven bubbles (Lancaster et al. 2021a,b).
For the Kroupa (2001) IMF, the total stellar mass formed for every SN progenitor star is m * ∼ 100 M , and the areal rate of SN explosions in quasi-steady state is Σ SFR /m * . For spherical momentum injection per SN of p * centered on the disk midplane, Ostriker & Shetty (2011) argued that the turbulent pressure ρv 2 z is expected to be comparable to the rate of vertical momentum flux injected on either side of the disk, P turb = (p * /4)(Σ SFR /m * ). Since p * is insensitive to the environment (both density and metallicity; e.g., Kim & Ostriker 2015a;Kim et al. 2017aKim et al. , 2023, the turbulent feedback yield is expected to be nearly constant (this is in stark contrast to the thermal feedback yield). The fit to the TIGRESS-classic suite gives (Ostriker & Kim 2022) log 10 (P turb,2p /k B ) = 0.96 log 10 (Σ SFR ) + 6.17 (24a) log 10 Υ turb = −0.06 log 10 (P DE /k B ) + 2.81.(24b) Figure 11(b) plots P turb vs. Σ SFR , showing the expected near linear relationship (Equation 24a). The turbulent feedback yield shown in Figure 11(e) is consistent with the shallow dependence on P DE seen in Equation 24b. We note that the current simulations have additional momentum injection by expanding H II regions as well as direct UV radiation pressure. Apparently, the contribution of UV in modulating global turbulent pressure is not significant. This strongly contrasts with the dominant role of UV in the destruction of molecular clouds (e.g., Kim et al. 2018Kim et al. , 2021.

Magnetic Stress
We find that the midplane magnetic stress and hence magnetic feedback yield (Figure 11(c) and (f)) is comparable to the turbulent kinetic component, for both models. Magnetic terms are determined by galactic dynamo action as a result of the interaction between turbulence (driven by feedback), galactic differential rotation, and buoyancy. The turbulent component of magnetic fields is directly related to the kinetic energy in turbulence and turbulent magnetic energy density quickly saturates at a level similar to kinetic energy density as long as the initial field is strong enough (Kim & Ostriker 2015b). Our initial field is purely azimuthal (along the y direction) and comparable to the final saturation level. Overall, the current simulations cover long-term evolution and result in a saturated state without a sign of further secular evolution in magnetic field strengths. 6

Total pressure and SFR prediction
Given the validity of vertical dynamical equilibrium and agreement of W with the simple weight estimator P DE (Equation 21), the PRFM theory postulates that the yield Υ tot = P tot /Σ SFR (calibrated from simulations) can be used to predict Σ SFR from P DE , which is calculated from large-scale galactic properties in observations. Summing up all pressure components, we obtain the total pressure support and the corresponding feedback yield. We find median Υ tot = 1500 km s −1 for model R8-4pc and 720 km s −1 for model LGR4-2pc, respectively.
As shown in Figure 12, the new simulation results are overall in good agreement with the TIGRESS-classic suite for the relation between Σ SFR and pressure or weight. In each panel, we directly compare our results with the fitting results from Ostriker & Kim (2022) for measured midplane pressure, measured weight, and es-6 The regular (mean) component of magnetic fields is maintained in our simulations as we include galactic differential rotation using the shearing box. In separate experiments without rotation or weak shear, we find much lower saturation level of magnetic fields and hence magnetic stress. We defer the detailed exploration and discussion of the magnetic field evolution to a separate work.
timated weight: log 10 (Σ SFR ) = 1.18 log 10 (P tot,2p /k B ) − 7.43 (25a) log 10 (Σ SFR ) = 1.17 log 10 (W/k B ) − 7.32 (25b) log 10 (Σ SFR ) = 1.21 log 10 (P DE,2p /k B ) − 7.66. (25c) We refrain from delivering a new fitting formula or making additional quantitative adjustments to the feedback yields given the limited parameter space covered in the present work. In the future, we will extend our parameter space survey, especially toward low metallicity regimes, to generalize the numerical calibration of feedback yield in the PRFM theory.

Summary of simulation results
We present first results from a new numerical framework that synthesizes the TIGRESS-classic computational model of the star-forming ISM (Kim & Ostriker 2017) with our recently developed non-equilibrium cooling and radiation (NCR) module (Kim et al. 2023). The detailed photochemical treatment and the effects of UV radiation from massive young stars are combined with the gravitational collapse/star formation and SN injection schemes implemented and tested in the TIGRESSclassic framework, in order to study the multiphase, turbulent, magnetized ISM self-consistently.
This paper considers two galactic conditions, one representing the solar neighborhood (R8) and the other a higher density/pressure environment (LGR4; close to the molecular gas weighted mean conditions in the PHANGS survey). We delineate the ISM properties, with a focus on the multiphase ISM distribution near the midplane (within a scale height). We then repeat the basic analysis done in Ostriker & Kim (2022) to test, validate, and calibrate the PRFM star formation theory. The key measured quantities from our analysis are summarized in Table 2, Table 6, and Figure 13.
We summarize the ISM phase distributions by mass and volume within |z| < 300 pc in Figure 13(a) and (b). Near the galactic midplane (within one scale height of the disk), the cold, unstable, and warm neutral medium (CNM, UNM, and WNM) occupies about 25%, 30%, and 35% by mass and 2%, 20%, and 50% by volume, respectively, in the solar neighborhood model R8. The warm ionized medium (WIM=WPIM+WCIM) contributes 8% and 10% by mass and volume, respectively, while the hot medium (Hot=WHIM+HIM) occupies about 15% of the volume with a negligible mass contribution. It is important to keep in mind that there are large amplitude temporal fluctuations (up to a factor 2) in these values, as indicated by the box (25-75 percentiles) and whisker (16-84 percentiles) in Figure 13. Moving to conditions of higher gas surface density, total pressure, and SFR with model LGR4, the mass contribution from the Cold (=CMM+CNM) component increases, while the volume filling factors remain more or less the same. For both models, the mass fractions of Cold increase with higher resolution at the expense of UNM and WNM, although the change in R8 is well within the time fluctuation level. The volume fractions are converged up to the temporal variation level. Figure 13(c) and (d) show the midplane pressure components and feedback yields for both models. The two different resolutions give converged results for both R8 and LGR4. The turbulent feedback yields are similar for R8 and LGR4, with a slightly decreasing trend toward higher pressure environment. The thermal feedback yield decreases as expected from R8 to LGR4, due to higher shielding of FUV radiation field in higher density environments. In an upcoming paper (N. Linzer et. al. in prep.), we will analyze the radiation field in depth to validate and calibrate the global attenuation model used in TIGRESS-classic (see Equation 22). The magnetic feedback yields in R8 are generally larger than those of LGR4; understanding the magnetic feedback yields requires further investigation of the galactic dynamo process, which in itself is a large and challenging area of research. As shown in Figure 12, the total feedback yields are quite similar to those reported in Ostriker & Kim (2022). For R8, Υ tot = 1500 km s −1 , and for LGR4, Υ tot = 720 km s −1 . Both models have similar σ eff,2p ≈ 12 − 13 km s −1 and σ z,turb,2p ∼ 7 − 8 km s −1 .
Finally, it is worth noting that the decrease in the WNM mass fraction from model R8 to model LGR4 is at least qualitatively consistent with theoretical expecta-tions (Ostriker & Kim 2022). The WNM mass fraction may be written as for f V,WNM the volume filling factor, where we have used P th = ρ w c 2 w for c w the warm-gas sound speed (which is insensitive to galactic environment), ρ w for the typical density in the warm medium near the midplane, and P tot ≡ ρ tot σ 2 eff for ρ tot the total midplane density. From Table 2 and Table 4, σ eff and f V,WNM are very similar between model R8 and LGR4, whereas Table 6 gives ∼ 30% lower ratio of thermal to total pressure for LGR4 than that for R8. About 30% decrease in the WNM mass fraction (Table 4 and Figure 13) is consistent with the expectation.

Comparison to Milky Way Empirical Constraints on ISM Phases
Our phase distribution is overall in good agreement with multiwavelength galactic observations. H I 21cm lines are the fundamental probe of the atomic ISM. An accurate determination of both gas column density and spin temperature requires H I absorption line measurements paired with emission lines. Generally speaking, WNM dominates 21 cm emission spectra, but WNM is extremely faint in absorption due to its low density and high spin temperature (which can be as high as the gas temperature due to radiative excitation by Lyα resonant scattering ;Wouthuysen 1952;Field 1959;Seon & Kim 2020). The detection of WNM (and UNM) in absorption requires highly sensitive absorption observations. There have been a number of sensitive absorption line surveys that determine mass fractions of different H I components (Heiles & Troland 2003;Roy et al. 2013;Murray et al. 2018). Using a simple radiative transfer model with multiple Gaussian components, a component de-tected in absorption with low or intermediate spin temperature (T s < 250 K or 250 K < T s < 1000 K) is considered to be CNM or UNM, while a component detected only in emission is WNM (with a small fraction detected in absorption with large spin temperature). The mass contribution to total H I column density of each component is roughly 30%, 20-30%, and 40-50% for CNM, UNM, and WNM, respectively, generally consistent among surveys. From Table 4, the mass fractions of the cold, unstable, and warm neutral medium in model R8 are ∼ 0.3, 0.3, 0.4, generally consistent with current empirical constraints for the solar neighborhood to the extent that they are available.
The observational measurement of thermal pressure using C I fine structure lines shows a lognormal distribution with a mean at P th,CNM ∼ 4000 K cm −3 and an rms dispersion of 0.175 dex (Jenkins & Tripp 2011). We obtain the mass-weighted pressure PDF in R8 with mean and standard deviation of ∼ 1500 K cm −3 and 0.27 dex for CNM, ∼ 3700 K cm −3 and 0.35 dex for UNM, and ∼ 4000 K cm −3 and 0.36 dex for WNM.
Observations of Hα and pulsar dispersion measures suggest that WIM forms a thick layer with scale height of ∼ 1 − 2 kpc (Reynolds 1989(Reynolds , 1991Taylor & Cordes 1993;Hill et al. 2008;Gaensler et al. 2008). One can deduce the volume-averaged midplane electron density n e ∼ 0.02−0.05 cm −3 (by using dispersion measures to pulsars with known distances) and filling factor f V,WIM ∼ 0.05 − 0.15 (by combining emission measures and dispersion measures) of the diffuse WIM (Kulkarni & Heiles 1987;Ferrière 2001;Gaensler et al. 2008). For the midplane number density of total gas n H ∼ 0.5 − 1 cm −3 (e.g., Boulares & Cox 1990;McKee et al. 2015), the mass fraction of WIM at the midplane is n e / n H < ∼ 10%. The mass fraction of WIM of f M,WIM ∼ 7% within |z| < 300 pc (see Table 4) is very much consistent with this empirical result.
Direct measurement of the hot gas in X-rays is difficult due to its low density. Also, significant diffuse soft Xray emission is from the Local Bubble (Cox & Reynolds 1987). Soft X-ray radiation from larger scales is presumably absorbed; for example, the band averaged absorption cross section at ∼ 0.25 keV is ∼ 10 −20 cm 2 H −1 (Snowden et al. 1990), yielding the mean free path ∼ 30 pc(n H /1 cm −3 ) −1 . Direct observational constraints on the larger scale distribution of likely pervasive hot gas in our Galaxy are still lacking.

Comparisons to Self-consistent Numerical Models of the Star-Forming ISM
Because the SFR and the ISM thermal and dynamical state co-regulate each other, one cannot be considered independently of the other. A theoretical model that explicitly addresses co-regulation, computing the SFR needed to maintain the thermal properties of the warm and cold ISM, was introduced by Ostriker et al. (2010); this and subsequent theoretical developments are summarized in Ostriker & Kim (2022).
Several groups have recently developed numerical frameworks that solve (magneto)hydrodynamics with cooling and heating, including stellar feedback (of various forms) from star clusters that are self-consistently formed via gravitational collapse.
Our own numerical studies began with a focus on just the warm and cold ISM, with feedback in the form of momentum injection and heating, both proportional to Σ SFR (Kim et al. 2011(Kim et al. , 2013Kim & Ostriker 2015b). These simulations, with a wide range of Σ gas , showed that a quasi-steady state is reached, validating vertical dynamical equilibrium. For a solar neighborhood model (QA10 in Kim et al. 2013), the values of Σ SFR ∼ 1.5 × 10 −3 M kpc −2 yr −1 and the midplane pressure (=weight) ∼ 10 4 k B cm −3 K were about a factor of two lower than those reported here (and from TIGRESS-classic) due to missing magnetic support and slightly weaker turbulence (H ∼ 80 pc vs. 220 pc and σ z,turb ∼ 5 km s −1 vs. 7 − 8 km s −1 ). Coincidentally, the total feedback yield (without magnetic contribution) in Kim et al. (2013) is similar to that of the current simulations as the fixed (p * /m * ) = 3 × 10 3 km s −1 adopted in the earlier work was higher than the effective (p * /m * ) eff ∼ 10 3 km s −1 realized via self-consistent expansion of SNe driven bubbles (Kim et al. 2017a). Kim & Ostriker (2017) introduced the TIGRESSclassic framework, with full treatment of the hot ISM. Direct comparison with the TIGRESS-classic suite results from Ostriker & Kim (2022) regarding SFR, pressures, and feedback yields show overall consistent results with the current work, modulo slightly larger value of Υ th and hence Υ tot here (see Section 4). The lack of local shielding of FUV in TIGRESS-classic tends to result in lower Cold mass fraction (f M,CNM+UNM ∼ 30% in TIGRESS-classic vs. f M,Cold ∼ f M,UNM ∼ 30% in TIGRESS-NCR). Inclusion of the ionizing radiation in TIGRESS-NCR converts significant WNM into WPIM (f M,WPIM ∼ 7−8%), which is similar to the value obtained from the post-processing of TIGRESS-classic (Kado-Fong et al. 2020).
Hennebelle & Iffrig (2014) and Colling et al. (2018) are similar to our earlier work (Kim et al. 2013;Kim & Ostriker 2015b) in terms of their SN feedback being mostly in the form of momentum injection without creating hot gas. The velocity dispersion in their models (which have Σ gas = 20 M pc −2 ) is about 5 − 7 km s −1 , which is slightly lower than both of our new models. The mag-netic fields tend to reduce SFR up to a factor of 2, with more reduction in the strong rotation case. Given that the magnetic feedback yield is about 30 − 40% of the total, Kim & Ostriker (2015b) found similarly higher SFR in non-magnetized cases. Comparisons between magnetized and unmagnetized cases using the TIGRESS-NCR framework will be investigated in a separate paper.
The SILCC framework first introduced by Walch et al. (2015) focuses on solar neighborhood ISM modeling, with particular emphasis on hydrodynamical evolution with a hydrogen and carbon chemistry network (Glover & Mac Low 2007a,b;, collectively called SGChem). Gatto et al. (2017) added sink particle treatments and star formation via gravitational collapse in the SILCC framework. They emphasized the role of stellar winds shutting off further accretion after sink formation. They found the resulting SFR surface density and ISM properties (mostly focused on hot gas filling factor) are sensitive functions of the density threshold for sink particle formation. The highest density threshold model (n thresh = 10 4 cm −3 ) yields Σ SFR ∼ 10 −3 M kpc −2 yr −1 , while the low density threshold model (n thresh = 10 2 cm −3 ) experiences a strong initial burst of star formation with Σ SFR > 10 −2 M kpc −2 yr −1 . Peters et al. (2017) included radiation transfer for ionizing UV (without radiation pressure and with a constant FUV background), with the same treatment of SNe and stellar winds, and a high density threshold. Their models with and without ionizing radiation (with both SNe and stellar winds) show similar Σ SFR ∼ 10 −3 M kpc −2 yr −1 but the inclusion of UV radiation gives smaller f V,h of ∼ 20 − 30%, larger warm gas filling factor, and reduced H 2 gas mass (about a factor two) at the end of their simulation (∼ 70 Myr). However, these quantities were still evolving, and the short runtime of their simulations makes it unclear whether the reported values are representative values in the statistical steady state of these models. The SFR obtained by Peters et al. (2017) in their simulations with SNe, stellar winds, and ionizing radiation is similar to what we obtain here for the R8 model, Σ SFR = 3 × 10 −3 M kpc −2 yr −1 .
Recently, Rathjen et al. (2021) conducted simulations using the SILCC framework with a more comprehensive feedback model including SNe, stellar winds, UV radiation, and CRs, as well as magnetic fields. By systematically turning on and off each feedback process, they found a progressive decrease in Σ SFR , f V,h , and cold gas mass fraction with more feedback. The impact of CRs is not significant (given the short evolution time of ∼ 100 Myr), and the model with SN, stellar winds, and radiation (called SWR) shows Σ SFR ∼ 1.5 − 2 × 10 −3 M kpc −2 yr −1 , similar to what we find and to observations. Within |z| < 250 pc, their SWR model shows f V,h ∼ 50% (35% with CRs) and f M,cold ∼ 50%; both are larger than what we find here. One potential reason is that their FUV radiation was assumed to be constant over time so that the thermal balance in the volume filling warm and cold ISM may not be fully self-consistent. EUV radiation was transferred using a tree-based backward ray tracing method (Wünsch et al. 2021), which is inherently less accurate than the direct ray tracing method we adopt here, especially behind regions of strong shielding (pervasive for EUV due to the huge cross section of neutral hydrogen against ionizing radiation). Finally, due to short evolution time (t < 100 Myr), their measurements include an initial burst period (25 Myr < t < 100 Myr), which may bias the hot gas filling factor toward higher values. Hu et al. (2021) developed a local simulation that handles time-dependent hydrogen chemistry on-the-fly using a chemistry network based on SGChem, and explored the effect of metallicity. Their radiation treatment is approximate: the (spatially-constant) unattenuated UV radiation field and CR ionization rate are scaled by recent star formation, with a local attenuation factor for FUV radiation applied using a tree-based method . Photoionization is treated using an iterative Strömgren sphere approach (Hu et al. 2017). Although properties of the ISM phase structure from this simulation were not explicitly discussed, Σ SFR ∼ 2 − 3 × 10 −3 M kpc −2 yr −1 and the mass fraction of warm ionized medium (∼ 5 − 10%) for the solar metallicity model are consistent with observations and our results.

Future Perspectives
The new simulation framework, TIGRESS-NCR, presented in this paper provides a promising tool for modeling the star-forming ISM. The main advance from the TIGRESS-classic framework is including direct UV radiation transfer and explicit chemical abundance calculations. These extensions allow us to examine more detailed aspects of ISM physics, and enable us to explore new parameter space beyond the conditions that apply in normal, low-redshift spiral galaxies like the Milky Way. One immediate application is to explore low metallicity environments that are common in local dwarfs and prevalent in all galaxies at high redshifts. Effects of metallicity on species abundances and the CO-to-H 2 conversion factor have been studied in previous work, with Gong et al. (2018Gong et al. ( , 2020 postprocessing the TIGRESS-classic suite with six-ray radiation transfer and steady-state chemistry, and Hu et al. (2021Hu et al. ( , 2022 using a tree-based shielding column calculation with time-dependent hydrogen chemistry combined with steady state carbon/oxygen chemistry. Given the more accurate methods for UV radiation transfer implemented in the TIGRESS-NCR framework, it will be very interesting to make comparisons with these works employing approximate radiation transfer. Also, the capability of modeling time-dependent H 2 chemistry will be an important tool in understanding observed chemical abundances (see Godard et al. 2023, for CH + and hence warm diffuse H 2 abundances), although higher resolution than is possible in the present simulations may be needed for many applications.
With a suite of simulations at varying metallicity, we can extend the theoretical understanding of SFR/ISM co-regulation to the low-metallicity regime, where the thermal feedback yield (and therefore thermal pressure) to become larger than other components because radiation easily propagates over large distances. This extension of the PRFM theory will be critical in developing a subgrid star formation recipe for large scale cosmological simulations. Applying the TIGRESS-NCR framework to study regions with strong spiral structures will be straightforward, since the TIGRESS-classic framework has already been successfully used for models of this kind (Kim et al. 2020b).
Although TIGRESS-NCR represents a significant advance in resolving and modeling key physical processes, there is still more to be done. First, we do not explicitly model CR transport. Currently, we only include ionization and heating by low-energy CRs, with the background value scaled with Σ SFR and Σ gas and attenuated in high density environments (see 2.3). This is a physically and empirically motivated prescription but lacks quantitative calibration from direct numerical modeling and ignores the dynamical effect of CRs. Full CR transport should include advection by the gas, streaming along magnetic field lines at the (ion) Alfven speed, and diffusion by scattering off of MHD waves that are likely self-generated for GeV and lower energies (Armillotta et al. 2021(Armillotta et al. , 2022. TIGRESS-NCR provides a unique laboratory for CR transport modeling as our framework produces a turbulent, multiphase ISM with realistic magnetic field and ionization structure as well as realistic, high-velocity hot galactic winds. Although ∼GeV CRs dominate the total energy budget and are expected to be dynamically important (Girichidis et al. 2018), low-energy CRs are responsible for ionization in most of the ISM's mass. Therefore, spectrally resolved CR transport is necessary (e.g., Girichidis et al. 2020Girichidis et al. , 2022. Thermal conduction, which is not included in our current framework, can alter the hot gas properties. The conductive heat transport from hot gas (created by SNe) to the warm/cold ISM leads to evaporation, although conductivity may be suppressed perpendicular to the magnetic field (Braginskii 1965). To the extent that it can act, conductive evaporation maintains the hot gas pressure while increasing its mass and decreasing its temperature (e.g., El-Badry et al. 2019). Conduction could certainly alter the observable properties of diffuse X-ray emission from the hot gas, and potentially change the hot gas mass fraction and volume filling factor.

ACKNOWLEDGMENTS
We are grateful to the referee for the timely and helpful report.