Introducing TIGRESS-NCR. I. Coregulation of the Multiphase Interstellar Medium and Star Formation Rates

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 et al. 2008; Stone & Gardiner 2009), employing a high-order Godunov method combined with a constrained transport algorithm (Gardiner & Stone 2008).

The conservation equations for gas mass, momentum, and total energy are, respectively,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{v}}\right)=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot \left[\rho {\boldsymbol{v}}{\boldsymbol{v}}+\left(P+{P}_{B}\right){\mathbb{I}}-\displaystyle \frac{{\boldsymbol{B}}{\boldsymbol{B}}}{4\pi }\right]\\ \quad =-\rho {\rm{\nabla }}{\rm{\Phi }}-2{\boldsymbol{\Omega }}\times (\rho {\boldsymbol{v}})+{{\boldsymbol{f}}}_{\mathrm{rad}},\end{array}\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {{ \mathcal E }}_{\mathrm{tot}}}{\partial t} & + & {\rm{\nabla }}\cdot \left[{\boldsymbol{v}}({{ \mathcal E }}_{\mathrm{tot}}+P+{P}_{B})-\displaystyle \frac{{\boldsymbol{B}}({\boldsymbol{B}}\cdot {\boldsymbol{v}})}{4\pi }\right]\\ & & =\,{ \mathcal G }-{ \mathcal L }-(\rho {\boldsymbol{v}})\cdot {\rm{\nabla }}{\rm{\Phi }}+{\boldsymbol{v}}\cdot {{\boldsymbol{f}}}_{\mathrm{rad}}.\end{array}\end{eqnarray} \tag{ 3 }$$

The magnetic field evolution is governed by the induction equation without explicit resistivity (ideal MHD):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times \left({\boldsymbol{v}}\times {\boldsymbol{B}}\right),\end{eqnarray} \tag{ 4 }$$

while B must obey the divergence-free constraint

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{B}}=0.\end{eqnarray} \tag{ 5 }$$

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 hydrogen 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; and ${{ \mathcal E }}_{\mathrm{tot}}=\rho {v}^{2}/2+P/(\gamma -1)+{P}_{B}$ is the total energy density, where γ = 5/3 is the adiabatic index.

We explicitly follow nonequilibrium abundances of molecular (${x}_{{{\rm{H}}}_{2}}$) and ionized (${x}_{{{\rm{H}}}^{+}}$) hydrogen by solving the transport of abundances with source terms,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{s}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{s}{\boldsymbol{v}})=\rho {{ \mathcal C }}_{s},\end{eqnarray} \tag{ 6 }$$

where ρ_s = ρ x_s is the mass density of species s (H₂ or H⁺), and ${{ \mathcal 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),

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{{\rm{\Phi }}}_{\mathrm{sg}}=4\pi G(\rho +{\rho }_{\mathrm{sp}}),\end{eqnarray} \tag{ 7 }$$

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 (${ \mathcal G }-{ \mathcal 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₀, i.e., ${\boldsymbol{\Omega }}={{\rm{\Omega }}}_{0}\hat{{\boldsymbol{z}}}$, and we assume the galactic rotation curve is flat, i.e., the shear parameter $q\equiv -d\mathrm{ln}{\rm{\Omega }}/d\mathrm{ln}R=1$. As a result, ${{\rm{\Phi }}}_{\mathrm{tidal}}(x)=-q{{\rm{\Omega }}}_{0}^{2}{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 TIGRESS-classic, 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 (${{ \mathcal E }}_{\mathrm{rad}};$ see Section 2.3), and the CR ionization rate (ξ_cr). Explicitly, we have

$$\begin{eqnarray}&&{{ \mathcal C }}_{s}\equiv {{ \mathcal C }}_{s}({n}_{{\rm{H}}},T,{x}_{s},{Z}_{d}^{{\prime} },{Z}_{g}^{{\prime} },{{ \mathcal E }}_{\mathrm{rad}},{\xi }_{\mathrm{cr}}),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{ \mathcal G }\equiv {n}_{{\rm{H}}}{\rm{\Gamma }}({n}_{{\rm{H}}},T,{x}_{s},{Z}_{d}^{{\prime} },{Z}_{g}^{{\prime} },{{ \mathcal E }}_{\mathrm{rad}},{\xi }_{\mathrm{cr}}),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{ \mathcal L }\equiv {n}_{{\rm{H}}}^{2}{\rm{\Lambda }}({n}_{{\rm{H}}},T,{x}_{s},{Z}_{d}^{{\prime} },{Z}_{g}^{{\prime} },{{ \mathcal E }}_{\mathrm{rad}},| {dv}/{dr}| ),\end{eqnarray} \tag{ 10 }$$

where ${Z}_{g}^{{\prime} }$ and ${Z}_{d}^{{\prime} }$ 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}^{{\prime} }={Z}_{d}^{{\prime} }=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}^{{\prime} }$ and ${Z}_{d}^{{\prime} }$, 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}^{{\prime} }$ and ${Z}_{d}^{{\prime} }$ except for gas with very low metal and dust content in the early universe. ⁶

As detailed in Kim et al. (2023), we note that the heating and cooling functions that we adopt follow Wolfire et al. (1995, 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 time dependent, they also make it possible to test the extent to which the predicted thermal equilibrium state is actually reached.

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 the 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³ 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 as follows: (i) a cell's density is higher than a threshold Larson–Penston density depending on sound speed and numerical resolution; (ii) the cell is at a local gravitational potential minimum; and (iii) flows along all three Cartesian directions converge toward the cell. Each particle represents a star cluster (with typical mass >10³ 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, 2014) to determine SN rates for each star cluster, assuming a Kroupa (2001) IMF and Geneva evolutionary tracks for nonrotating 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 3Δx. If the calculated gas mass is less than the shell formation mass ${M}_{\mathrm{sf}}=1540\,{M}_{\odot }{({n}_{\mathrm{amb}}/{\mathrm{cm}}^{-3})}^{-0.33}$ when a remnant becomes radiative (Kim & Ostriker 2015a), a total of ${E}_{\mathrm{SN}}={10}^{51}\,\mathrm{erg}$ energy is injected with a thermal-to-kinetic energy ratio consistent with that of the Sedov–Taylor stage of evolution, after averaging out the feedback injection region. If the Sedov–Taylor stage is unresolved (i.e., M_SNR > M_sf), a total of ${p}_{\mathrm{SNR}}=2.8\times {10}^{5}\,{M}_{\odot }\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{({n}_{\mathrm{amb}}/{\mathrm{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.

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 the 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₂ and CO, and ionize C. To compute the dissociating radiation field for H₂, we apply the Draine & Bertoldi (1996) self-shielding function to the LW band, using the H₂ column density integrated from each source to each cell. The LyC photons (also referred to as extreme-UV, EUV) ionize neutral hydrogen (H and H₂) 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,

$$\begin{eqnarray}&&\hat{{\boldsymbol{k}}}\cdot {\rm{\nabla }}{I}_{j}=-{\alpha }_{j}{I}_{j},\end{eqnarray} \tag{ 11 }$$

for each frequency bin j along a set of rays. Here, I_j is the intensity, $\hat{{\boldsymbol{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⁴ 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 four 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, ${{{ \mathcal E }}_{\mathrm{rad}}}_{,j}$, and flux, ${{{\boldsymbol{F}}}_{\mathrm{rad}}}_{,j}$, in each cell are obtained by summing over contributions from all rays passing through it. We then have ${{\boldsymbol{f}}}_{\mathrm{rad}}={\sum }_{j}\rho {\kappa }_{j}{{{\boldsymbol{F}}}_{\mathrm{rad}}}_{,j}/c$, which we use to update the gas momentum and corresponding kinetic energy density; the values of ${{{ \mathcal E }}_{\mathrm{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 ($\hat{{\boldsymbol{x}}}$) edges of the computational domain are offset by the distance the boundaries have sheared in the local-azimuthal ($\hat{{\boldsymbol{y}}}$) 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 z-boundary of the computational domain or the optical depth measured from the source is larger than ${\tau }_{\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 time step is prohibitive. To reduce the cost of ray tracing without losing accuracy too much, we adopt three modifications: (i) we perform ray tracing at intervals Δt_2p based on the Courant–Friedrichs–Lewy time step for the cold and warm gas at T < 2.0 × 10⁴ K, or whenever a new radiation source is created, or whenever an existing radiation source is removed; (ii) we put a hard limit on the maximum horizontal distance a ray may propagate from each source, which we denote as ${d}_{\mathrm{xy},\max }$; (iii) 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 (ii) reduces the cost to an acceptable level without significantly affecting the radiation field solution in the midplane region, provided ${d}_{\mathrm{xy},\max }$ is large enough (see Appendix A.2). The condition (iii) 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}_{\mathrm{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 Appendix 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 \theta =\hat{{\boldsymbol{k}}}\cdot \hat{{\boldsymbol{z}}}$. 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 (iii), 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 ${\xi }_{\mathrm{cr},0}=2\times {10}^{-16}\,{{\rm{s}}}^{-1}{{\rm{\Sigma }}}_{\mathrm{SFR},40}^{{\prime} }/{{\rm{\Sigma }}}_{\mathrm{gas}}^{{\prime} }$, where ${{\rm{\Sigma }}}_{\mathrm{SFR},40}^{{\prime} }$ and ${{\rm{\Sigma }}}_{\mathrm{gas}}^{{\prime} }$ are the SFR surface density measured from stars formed in the last 40 Myr and instantaneous gas surface density relative to solar-neighborhood values Σ_SFR = 2.5 × 10⁻³ M_⊙ kpc⁻² yr⁻¹ and Σ_gas = 10 M_⊙ pc⁻². The local CR ionization rate is then set to

$$\begin{eqnarray}{\xi }_{\mathrm{cr}}=\left\{\begin{array}{ll}{\xi }_{\mathrm{cr},0} & {\rm{if}}\,{N}_{\mathrm{eff}}\leqslant {N}_{0};\\ {\xi }_{\mathrm{cr},0}{\left(\displaystyle \frac{{N}_{\mathrm{eff}}}{{N}_{0}}\right)}^{-1} & {\rm{if}}\,{N}_{\mathrm{eff}}\gt {N}_{0},\end{array}\right.\end{eqnarray} \tag{ 12 }$$

where N_eff is the effective shielding column density and N₀ = 9.35 × 10²⁰ cm⁻². 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}_{\mathrm{eff}}={10}^{21}\,{\mathrm{cm}}^{-2}{Z}_{d}{{\prime} }^{-1}\mathrm{ln}({{ \mathcal E }}_{\mathrm{PE},\mathrm{unatt}}/{{ \mathcal E }}_{\mathrm{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/{{\rm{\Sigma }}}_{\mathrm{gas}}^{{\prime} }$ in setting ξ_cr,0 represents attenuation under average conditions and is motivated by Wolfire et al. (2003). The additional attenuation at high column densities in Equation (12) is motivated by observations of column-density-dependent CR ionization rate (Neufeld & Wolfire 2017).

In the MHD integrator, we transport molecular and ionized hydrogen using passive scalars (${x}_{{{\rm{H}}}_{2}}$ and ${x}_{{{\rm{H}}}^{+}}$, respectively) without source terms as implemented in Athena. We obtain the atomic hydrogen abundance from the closure relation ${x}_{{\rm{H}}}=1-2{x}_{{{\rm{H}}}_{2}}-{x}_{{{\rm{H}}}^{+}}$. We then update the temperature and abundances in an operator-split manner by solving two ordinary differential equations (ODEs):

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=K\displaystyle \frac{{{dT}}_{\mu }}{{dt}}={ \mathcal G }-{ \mathcal L },\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dx}}_{s}}{{dt}}={{ \mathcal C }}_{s},\,\,\,\,({\rm{s}}:\,{{\rm{H}}}_{2}\,\mathrm{and}\,{{\rm{H}}}^{+}),\end{eqnarray} \tag{ 14 }$$

where e = P/(γ − 1) is the internal energy density, T_μ ≡ T/μ for $\mu ={\mu }_{{\rm{H}}}/(1.1+{x}_{{\rm{e}}}-{x}_{{{\rm{H}}}_{2}})$ 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 timescales, 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 timescales.

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, ${ \mathcal G }-{ \mathcal L }$, with respect to temperature, which depends on the previous step's abundances and other information (see parameter dependence in Equations (9) and (10) as well as Section 4 of Kim et al. 2023). We then evaluate ${{ \mathcal 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⁴ K, or from helium and metals assuming CIE at T > 2 × 10⁴ 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.

• 1.  
  Molecular hydrogen: formation by grain catalysis; destruction by CR ionization, photodissociation, photoionization, and collisional dissociation.
• 2.  
  Ionized hydrogen: formation by photoionization, CR ionization, and collisional ionization of atomic hydrogen; destruction by radiative recombination and grain-assisted recombination.
• 3.  
  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.
• 4.  
  Heating: photoelectric effect on small grains by FUV photons; CR ionization of H and H₂; photoionization of H and H₂; formation, photodissociation, and UV pumping of H₂.
• 5.  
  Cooling:

  • (a)  
    Λ_hyd: collisionally excited Lyα resonance line from H; collisional ionization of H; collisional dissociation of H₂; rovibrational lines from H₂; bremsstrahlung and radiative/grain-assisted recombination of free electrons with H⁺.
  • (b)  
    Λ_others(T < 2 × 10⁴ K): fine-structure lines from C⁺, C, and O; rotational lines of CO; combined nebular lines in ionized gas (Λ_neb); grain-assisted recombination.
  • (c)  
    Λ_CIE(T > 3.5 × 10⁴ K): ion-by-ion CIE cooling table for He and metals from Gnat & Ferland (2012). Metal cooling is scaled linearly with ${Z}_{g}^{{\prime} }$.

  We smoothly transition from Λ_others to Λ_CIE at 2 × 10⁴ K < T < 3.5 × 10⁴ 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⁻⁵.

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 area-weighted 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)² × 6144 pc for R8, and (512 pc)² × 3072 pc for LGR4. We use ${d}_{\mathrm{xy},\max }=2048\,\mathrm{pc}$ and ε_PP = 0 for R8 and ${d}_{\mathrm{xy},\max }=1024\,\mathrm{pc}$ and ε_PP = 10⁻⁸ for LGR4. With these choices of numerical parameters for ray tracing, we found good convergence for the EUV radiation field everywhere, the FUV field near the midplane, and overall simulation outcomes (see Appendix A.2). We note that the selection of the optimal values of ${d}_{\mathrm{xy},\max }$ and ε_PP depends on the system's conditions, especially on ${Z}_{d}^{{\prime} }$ (which determines dust absorption).

We run each simulation with two different resolutions: 8 and 4 pc for R8, and 4 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 amplitudes of 10 and 15 km s⁻¹ for R8 and LGR4, respectively, along with randomly distributed initial star clusters that provide nonzero initial heating and SNe. The initial magnetic field is set to be azimuthal ${\boldsymbol{B}}={B}_{0}\hat{{\boldsymbol{y}}}$ with uniform plasma beta ${\beta }_{0}\equiv {P}_{\mathrm{th}}/({B}_{0}^{2}/8\pi )=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.

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,

$$\begin{eqnarray}&&{{{\rm{\Sigma }}}_{\mathrm{SFR}}}_{,10}\equiv \displaystyle \frac{{\sum }_{i}{M}_{* ,i}({t}_{\mathrm{age}}\lt {\rm{\Delta }}t)}{{L}_{x}{L}_{y}{\rm{\Delta }}t},\end{eqnarray} \tag{ 15 }$$

where M_*,i (t_age < Δt) is the total mass of star particles with age younger than Δt; (c) the effective vertical velocity dispersion,

$$\begin{eqnarray}&&{\sigma }_{\mathrm{eff}}\equiv {\left[\displaystyle \frac{\int \left(\rho {v}_{z}^{2}+P+{B}^{2}/(8\pi )-{B}_{z}^{2}/(4\pi )\right){dV}}{\int \rho {dV}}\right]}^{1/2};\end{eqnarray} \tag{ 16 }$$

and (d) the mass-weighted scale height of the gas,

$$\begin{eqnarray}&&H\equiv {\left(\displaystyle \frac{\int \rho {z}^{2}{dV}}{\int \rho {dV}}\right)}^{1/2}.\end{eqnarray} \tag{ 17 }$$

Here, the values of σ_eff and H are calculated only for the warm and cold gas with temperature T < 3.5 × 10⁴ 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., ${\sigma }_{z,\mathrm{turb}}\equiv {\left(\int \rho {v}_{z}^{2}{dV}/\int \rho {dV}\right)}^{1/2}$, and report this separately in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Global Properties at Saturation

Model	Σgas	${{{\rm{\Sigma }}}_{\mathrm{SFR}}}_{,10}$	H	σeff	σz,turb	tver	tdep
	(M⊙ pc−2)	(10−3 M⊙ kpc−2 yr−1)	(pc)	(km s−1)	(km s−1)	(Myr)	(Gyr)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
R8-4pc	${10.6}_{-0.2}^{+0.2}$	${2.8}_{-1.0}^{+1.5}$	${199.3}_{-44.3}^{+36.2}$	${12.0}_{-0.5}^{+1.0}$	${7.7}_{-0.7}^{+1.0}$	${23.3}_{-3.8}^{+8.0}$	${3.6}_{-1.0}^{+1.0}$
R8-8pc	${10.3}_{-0.2}^{+0.4}$	${2.8}_{-2.0}^{+3.5}$	${233.5}_{-57.9}^{+147.1}$	${13.5}_{-1.3}^{+2.9}$	${9.4}_{-2.1}^{+3.3}$	${24.3}_{-5.2}^{+16.7}$	${3.2}_{-1.4}^{+2.5}$
LGR4-2pc	${37.9}_{-0.9}^{+1.3}$	${34.8}_{-10.7}^{+10.4}$	${164.5}_{-47.1}^{+31.4}$	${13.4}_{-0.7}^{+0.7}$	${8.3}_{-0.8}^{+0.9}$	${17.8}_{-3.9}^{+6.7}$	${1.2}_{-0.2}^{+0.2}$
LGR4-4pc	${36.2}_{-0.6}^{+1.4}$	${29.0}_{-8.5}^{+20.3}$	${176.2}_{-66.2}^{+64.0}$	${14.6}_{-1.1}^{+1.2}$	${9.9}_{-1.2}^{+1.5}$	${15.4}_{-4.0}^{+9.6}$	${1.0}_{-0.1}^{+0.4}$

Note. Each column provides median values as well as the 16th to 84th percentile range, over t = 250–450 Myr for R8 and t = 250–350 Myr for LGR4. Column (2): gas surface density. Column (3): SFR surface density. Column (4): mass-weighted gas scale height. Column (5): effective vertical velocity dispersion. Column (6): turbulent component of vertical velocity dispersion. Column (7): vertical dynamical timescale. Column (8): gas depletion time. Note that Columns (4)–(7) are calculated for the warm and cold gas with temperature T < 3.5 × 10⁴ K. See text for definitions.

Download table as:  ASCIITypeset image

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 mesh-refined, restarted models are run for longer than one orbit time (or three to four star formation and feedback cycles) to obtain a fair statistical sample of states at higher resolution. In Table 2, we summarize means and standard deviations of the 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. ⁷

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 $\sim {(4\pi G{\rho }_{\mathrm{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. The quasi-periodic 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 and 300 Myr in R8-8pc).

Figure 2(a) shows the total energy injection rate per unit area by UV radiation ⁸ ${\dot{S}}_{\mathrm{UV}}\equiv {L}_{\mathrm{UV},\mathrm{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.

Zoom In Zoom Out Reset image size

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 ${\dot{S}}_{\mathrm{heat}}\equiv \int { \mathcal G }{dV}/({L}_{x}{L}_{y})$ shown in Figure 2(b) is the sum of hydrogen photoionization heating by LyC radiation (${\dot{S}}_{\mathrm{heat},\mathrm{PI}}/{\dot{S}}_{\mathrm{heat}}\sim 75 \%$), photoelectric heating from FUV (PE+LW) radiation on small dust grains (${\dot{S}}_{\mathrm{heat},\mathrm{PE}}/{\dot{S}}_{\mathrm{heat}}\sim 20 \%$), and CR ionization heating (${\dot{S}}_{\mathrm{heat},\mathrm{CR}}/{\dot{S}}_{\mathrm{heat}}\sim 1 \% \mbox{--}2 \% )$, plus a tiny contribution from H₂ heating (<0.1%). The global heating efficiency, defined as the ratio of the total heating rate to the UV radiation injection rate, is ${\dot{S}}_{\mathrm{heat},\mathrm{PI}+\mathrm{PE}}/{\dot{S}}_{\mathrm{UV}}\sim 5 \% \mbox{--}6 \%$, with individual efficiencies ${\dot{S}}_{\mathrm{heat},\mathrm{PE}}/{\dot{S}}_{\mathrm{FUV}}\sim 2 \%$ and ${\dot{S}}_{\mathrm{heat},\mathrm{PI}}/{\dot{S}}_{\mathrm{EUV}}\sim 15 \%$.

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 (${\dot{S}}_{\mathrm{SN}};$ Figure 2(d)). ${\dot{S}}_{\mathrm{SN}}$ is about two orders of magnitude smaller than the UV radiation injected (${\dot{S}}_{\mathrm{UV}};$ Figure 2(a)), and a factor ∼3 lower than the radiative heating rate ${\dot{S}}_{\mathrm{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 ${\dot{S}}_{\mathrm{SN}}$) accounts for the small excess of SN injection energy over the net cooling within the domain.

The instantaneous spatial distribution of the ISM's 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 Figures 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.

Zoom In Zoom Out Reset image size

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 emission measure (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 hydrogen, helium, and metals. In contrast, x_e ≈ 1 (dark green) in the warm ionized gas as electrons are mostly from photoionized hydrogen (and 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%; white-blue 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 the 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 ${ \mathcal 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 ${ \mathcal G }\mbox{--}{ \mathcal L }$ panel). The net cooling rate is highest at hot bubble boundaries (CIE cooling at T ∼ 10⁵ K). These interfaces where hot gas mixes with denser gas to become strongly radiative are important in bubble energetics.

As the 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-of-midplane 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.

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 nine 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 five 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 "warm-collisionally ionized" medium with a temperature cut at T = 1.5 × 10⁴ K, above which hydrogen can be collisionally ionized (see red dashed line in Figure 4(c)). We assign every cell to one of the nine phases exclusively.

Zoom In Zoom Out Reset image size

Table 3. Phase Definition

Name	Temperature	Abundance	Shorthand
Cold molecular medium a	T < 6 × 103 K	${x}_{{{\rm{H}}}_{2}}\gt 0.25$	CMM
Cold neutral medium b	T < 500 K	xH > 0.5	CNM
Unstable neutral medium	500 K < T < 6 × 103 K	xH > 0.5	UNM
Unstable ionized medium c	T < 6 × 103 K	${x}_{{{\rm{H}}}^{+}}\gt 0.5$	${\mathtt{UIM}}$
Warm neutral medium	6 × 103 K < T < 3.5 × 104 K	xH > 0.5	WNM
Warm photoionized medium	6 × 103 K < T < 1.5 × 104 K	${x}_{{{\rm{H}}}^{+}}\gt 0.5$	${\mathtt{WPIM}}$
Warm collisionally ionized medium	1.5 × 104 K < T < 3.5 × 104 K	${x}_{{{\rm{H}}}^{+}}\gt 0.5$	${\mathtt{WCIM}}$
Warm-hot ionized medium	3.5 × 104 K < T < 5 × 105 K	⋯	${\mathtt{WHIM}}$
Hot ionized medium	5 × 105 K < T	⋯	${\mathtt{HIM}}$

Notes.

^a This includes unstable temperature range but is dominated by cold. ^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.

Download table as:  ASCIITypeset image

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. 2023, in preparation).

Table 4. Mass Fractions and Volume Filling Factors with ∣z∣ < 300 pc

Model	Cold	UNM	WNM	WPIM	WCIM	WHIM	HIM
Mass fractions
R8-4pc	${27.4}_{-10.3}^{+4.8}$	${28.8}_{-3.1}^{+3.4}$	${34.9}_{-5.9}^{+10.6}$	${7.5}_{-3.1}^{+5.6}$	${0.2}_{-0.1}^{+0.2}$	${0.1}_{-0.0}^{+0.1}$	${0.05}_{-0.02}^{+0.05}$
R8-8pc	${22.0}_{-10.8}^{+11.6}$	${29.1}_{-8.7}^{+5.6}$	${36.0}_{-9.6}^{+15.0}$	${6.7}_{-5.4}^{+10.6}$	${0.3}_{-0.1}^{+0.3}$	${0.2}_{-0.1}^{+0.2}$	${0.07}_{-0.03}^{+0.09}$
LGR4-2pc	${37.3}_{-5.3}^{+3.2}$	${27.4}_{-0.9}^{+1.3}$	${27.2}_{-3.3}^{+3.1}$	${7.3}_{-1.7}^{+2.3}$	${0.2}_{-0.1}^{+0.1}$	${0.1}_{-0.0}^{+0.1}$	${0.07}_{-0.02}^{+0.03}$
LGR4-4pc	${31.6}_{-4.9}^{+3.8}$	${30.0}_{-1.3}^{+1.5}$	${29.9}_{-3.5}^{+4.5}$	${7.3}_{-2.1}^{+3.6}$	${0.3}_{-0.1}^{+0.1}$	${0.2}_{-0.1}^{+0.1}$	${0.10}_{-0.03}^{+0.03}$
Volume-filling factors
R8-4pc	${1.6}_{-1.0}^{+0.5}$	${17.0}_{-6.0}^{+4.2}$	${48.2}_{-5.9}^{+6.9}$	${10.9}_{-4.4}^{+7.6}$	${1.4}_{-0.3}^{+0.3}$	${3.6}_{-0.9}^{+1.1}$	${12.5}_{-3.8}^{+15.5}$
R8-8pc	${1.5}_{-1.1}^{+1.9}$	${14.4}_{-9.3}^{+8.5}$	${44.5}_{-11.9}^{+11.8}$	${10.0}_{-7.5}^{+12.5}$	${1.7}_{-0.5}^{+0.7}$	${4.0}_{-1.4}^{+1.4}$	${17.5}_{-9.3}^{+20.7}$
LGR4-2pc	${2.5}_{-0.6}^{+0.7}$	${14.5}_{-2.0}^{+4.2}$	${51.3}_{-8.4}^{+10.4}$	${9.3}_{-2.1}^{+3.2}$	${1.1}_{-0.3}^{+0.2}$	${3.1}_{-1.0}^{+0.7}$	${15.4}_{-8.9}^{+6.6}$
LGR4-4pc	${2.1}_{-0.6}^{+1.0}$	${12.6}_{-1.4}^{+3.9}$	${50.9}_{-8.2}^{+6.1}$	${9.8}_{-3.3}^{+4.1}$	${1.8}_{-0.3}^{+0.3}$	${3.7}_{-0.7}^{+0.7}$	${18.3}_{-7.2}^{+4.9}$

Note. Each column shows the median and 16th and 84th percentile range over t = 250–450 Myr and t = 250–350 Myr for R8 and LGR4, respectively.

Download table as:  ASCIITypeset image

Separating the warm and cold gas into Cold, UNM+WNM, 2p, and WIM components, Table 5 shows the effective vertical velocity dispersion as defined by Equation (16) and the mass-weighted turbulent velocity dispersion only considering the ${P}_{\mathrm{turb}}=\rho {v}_{z}^{2}$ term for each component. Given that the neutral medium (especially, warmer component UNM+WNM) dominates the mass fraction (Table 4), σ_{eff,2 p} agrees well with the effective velocity dispersion of all warm and cold gas at T < 3.5 × 10⁴ 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 mass-dominating component. On the other hand, it shows that WIM tracers will typically overestimate the mass-weighted 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.

Table 5. Mass-weighted Vertical Velocity Dispersion

Model	σeff				σz,turb
	Cold	UNM+WNM	2p	WIM	Cold	UNM+WNM	2p	WIM
R8-4pc	${5.4}_{-0.9}^{+1.2}$	${12.6}_{-0.7}^{+0.6}$	${11.2}_{-0.6}^{+0.6}$	${18.4}_{-1.9}^{+2.1}$	${4.2}_{-0.6}^{+1.5}$	${7.6}_{-0.7}^{+1.1}$	${6.9}_{-0.7}^{+1.0}$	${13.0}_{-2.0}^{+2.5}$
R8-8pc	${5.8}_{-0.7}^{+1.8}$	${13.1}_{-1.1}^{+2.2}$	${12.2}_{-1.5}^{+2.2}$	${21.1}_{-3.4}^{+10.6}$	${4.8}_{-1.1}^{+1.5}$	${8.6}_{-1.8}^{+2.1}$	${8.1}_{-1.8}^{+2.5}$	${16.2}_{-4.1}^{+12.3}$
LGR4-2pc	${6.1}_{-0.7}^{+1.1}$	${14.8}_{-1.2}^{+0.6}$	${12.9}_{-0.7}^{+0.7}$	${18.7}_{-0.8}^{+1.1}$	${5.3}_{-0.7}^{+1.4}$	${8.4}_{-0.7}^{+1.1}$	${7.7}_{-0.7}^{+1.1}$	${12.9}_{-1.3}^{+1.6}$
LGR4-4pc	${7.7}_{-1.2}^{+3.3}$	${15.9}_{-1.4}^{+1.4}$	${13.8}_{-1.2}^{+1.6}$	${21.0}_{-1.4}^{+2.0}$	${6.7}_{-1.8}^{+3.2}$	${10.1}_{-1.3}^{+1.8}$	${9.2}_{-1.1}^{+2.1}$	${14.8}_{-1.6}^{+2.5}$

Note. Each column shows the median and 16th and 84th percentile range over t = 250–450 Myr for R8 and t = 250–350 Myr for LGR4.

Download table as:  ASCIITypeset image

Figures 5 and 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 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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}_{{{\rm{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}_{{\rm{e}}}-0.5{x}_{{{\rm{H}}}_{2}}){n}_{{\rm{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.

Zoom In Zoom Out Reset image size

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 spatially nonuniform radiation (see Equations (9) and (10)). 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 cooling-dominated and heating-dominated regions. Given ξ_cr = 2.9 × 10⁻¹⁶ s⁻¹ for the background gas and the median value of χ_FUV = 0.87, in Figure 7(b) we plot an equilibrium curve as a 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³ cm⁻³ 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⁻³ and P/k_B ∼ 5 × 10³ cm⁻³ K and the minimum thermal equilibrium pressure of the CNM (the bulk net cooling region shown as blue) at n_H ∼ 5 cm⁻³ and P/k_B ∼ 1 × 10³ cm⁻³ K. Although the ionized gas has a very wide range of thermal pressure, the mean is P/k_B = 7.9 × 10³ cm⁻³ 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 adiabatic 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. However, 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.

Zoom In Zoom Out Reset image size

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 (${ \mathcal W }$):

$$\begin{eqnarray}&&{\rm{\Delta }}\unicode{x03008}{P}_{\mathrm{tot}}\unicode{x03009}\equiv {\rm{\Delta }}\unicode{x03008}{P}_{\mathrm{th}}+{P}_{\mathrm{turb}}+{{\rm{\Pi }}}_{\mathrm{mag}}\unicode{x03009}+{\rm{\Delta }}{P}_{\mathrm{rad}}={ \mathcal W },\end{eqnarray} \tag{ 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}_{\mathrm{turb}}\equiv \rho {v}_{z}^{2}$, and magnetic stress (magnetic pressure + tension) ${{\rm{\Pi }}}_{\mathrm{mag}}\equiv ({B}_{x}^{2}+{B}_{y}^{2}-{B}_{z}^{2})/(8\pi )$. Note that the mean or turbulent magnetic stress (${{\rm{\Pi }}}_{\overline{B}}$ or Π_{δ B} ) is, respectively, defined by substituting for B with the mean $\overline{{\boldsymbol{B}}}\equiv \unicode{x03008}{\boldsymbol{B}}\unicode{x03009}$ or turbulent $\delta {\boldsymbol{B}}\equiv {\boldsymbol{B}}-\overline{{\boldsymbol{B}}}$ component. The radiation pressure and weight terms can be defined toward either the upper or lower disk, ${\rm{\Delta }}{P}_{\mathrm{rad}}={\int }_{0}^{\pm {L}_{z}/2}\unicode{x03008}{f}_{\mathrm{rad},{\rm{z}}}\unicode{x03009}{dz}$ and ${ \mathcal W }={\int }_{0}^{\pm {L}_{z}/2}\unicode{x03008}\rho \partial {\rm{\Phi }}/\partial z\unicode{x03009}{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 ${{ \mathcal W }}_{\mathrm{ext}}$ (due to stellar disk and dark matter halo), and gas weight from self-gravity ${{ \mathcal W }}_{\mathrm{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

$$\begin{eqnarray}&&{\unicode{x03008}q\unicode{x03009}}_{\mathrm{ph}}\equiv \displaystyle \frac{\int \int q{\rm{\Theta }}(\mathrm{ph}){dxdy}}{{L}_{x}{L}_{y}},\end{eqnarray} \tag{ 19 }$$

where Θ(ph) = 1 if the temperature and abundance conditions in Table 3 are 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 $\unicode{x03008}{P}_{\mathrm{tot}}\unicode{x03009}={\sum }_{\mathrm{ph}}{\unicode{x03008}{P}_{\mathrm{tot}}\unicode{x03009}}_{\mathrm{ph}}$. Individual pressure components ($\unicode{x03008}{P}_{\mathrm{th}}\unicode{x03009}$, 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 by ${\tilde{P}}_{\mathrm{tot},\mathrm{ph}}\equiv {\unicode{x03008}{P}_{\mathrm{tot}}\unicode{x03009}}_{\mathrm{ph}}/{f}_{{\rm{A}},\mathrm{ph}}$ (and similarly for ${\tilde{P}}_{\mathrm{th},\mathrm{ph}}$, etc.). We can then write $\unicode{x03008}{P}_{\mathrm{tot}}\unicode{x03009}={\sum }_{\mathrm{ph}}{f}_{{\rm{A}},\mathrm{ph}}{\tilde{P}}_{\mathrm{tot},\mathrm{ph}}$. If the typical values of the total pressure for 2p, WIM, and Hot are comparable with each other, we have $\unicode{x03008}{P}_{\mathrm{tot}}\unicode{x03009}\approx {\tilde{P}}_{\mathrm{tot},{\rm{X}}}{\sum }_{\mathrm{ph}}{f}_{{\rm{A}},\mathrm{ph}}={\tilde{P}}_{\mathrm{tot},{\rm{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

$$\begin{eqnarray}&&{\rm{\Delta }}\unicode{x03008}{P}_{\mathrm{tot}}\unicode{x03009}\approx {\tilde{P}}_{\mathrm{tot},2{\rm{p}}}={\tilde{P}}_{\mathrm{th},2{\rm{p}}}+{\tilde{P}}_{\mathrm{turb},2{\rm{p}}}+{\tilde{{\rm{\Pi }}}}_{\mathrm{mag},2{\rm{p}}}={ \mathcal W }.\end{eqnarray} \tag{ 20 }$$

We note that weight (and radiation pressure) is vertically integrated over all phases, and that the (time-averaged) 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 the 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 the 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 (${ \mathcal W }$) and the commonly used weight estimator (which assumes that the stellar disk is thicker than the gaseous disk; Ostriker & Kim 2022):

$$\begin{eqnarray}&&{P}_{\mathrm{DE},2{\rm{p}}}\equiv \displaystyle \frac{\pi G{{\rm{\Sigma }}}_{\mathrm{gas}}^{2}}{2}+{{\rm{\Sigma }}}_{\mathrm{gas}}{\left(2G{\rho }_{\mathrm{sd}}\right)}^{1/2}{\sigma }_{\mathrm{eff},2{\rm{p}}},\end{eqnarray} \tag{ 21 }$$

constructed from observables (e.g., Sun et al. 2020). We note that σ_{eff,2 p} , presented in Table 5, is the mass-weighted mean for the 2p phase over the entire domain (not the midplane measure). This kinetic (thermal and turbulent) velocity dispersion is a direct observable given line emission from the neutral (atomic and molecular) gas, and then σ_{eff,2 p} 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).

Zoom In Zoom Out Reset image size

Figure 10 plots P_tot,2p as a function of (a) P_tot,hot, (b) ${ \mathcal 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 pressure. Thermal pressure dominates in WIM and Hot, while thermal pressure is the smallest component in 2p, with magnetic and turbulent components being comparable. ${P}_{\mathrm{tot},2{\rm{p}}}\approx { \mathcal 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.

Zoom In Zoom Out Reset image size

Table 6. Midplane Pressure and Weight

Model	Pth,2p	Pturb,2p	Πmag,2p	Ptot,2p	${P}_{\mathrm{tot},\mathrm{WIM}}$	Ptot,hot	Ptot	${ \mathcal W }$	PDE,2p
R8-4pc	${4.6}_{-0.8}^{+1.5}$	${7.4}_{-1.9}^{+2.5}$	${7.7}_{-2.7}^{+2.4}$	${20.2}_{-2.5}^{+2.9}$	${19.4}_{-3.6}^{+3.7}$	${15.1}_{-3.8}^{+5.2}$	${18.5}_{-3.4}^{+3.1}$	${17.6}_{-1.9}^{+2.8}$	${20.1}_{-1.5}^{+1.2}$
R8-8pc	${4.3}_{-1.2}^{+3.5}$	${7.0}_{-2.0}^{+5.1}$	${8.5}_{-2.5}^{+3.0}$	${20.6}_{-4.8}^{+6.7}$	${20.5}_{-5.4}^{+6.1}$	${17.8}_{-6.3}^{+6.2}$	${20.5}_{-5.7}^{+5.7}$	${19.1}_{-2.0}^{+5.6}$	${21.5}_{-2.4}^{+2.7}$
LGR4-2pc	${2.1}_{-0.1}^{+0.4}$	${5.6}_{-1.3}^{+3.6}$	${4.4}_{-2.2}^{+1.4}$	${13.0}_{-2.6}^{+1.5}$	${9.1}_{-1.9}^{+2.5}$	${9.5}_{-2.1}^{+4.5}$	${10.5}_{-2.2}^{+2.8}$	${10.7}_{-1.4}^{+1.1}$	${10.0}_{-0.4}^{+0.3}$
LGR4-4pc	${2.1}_{-0.7}^{+0.9}$	${6.2}_{-2.9}^{+1.5}$	${5.2}_{-3.1}^{+2.9}$	${12.4}_{-2.4}^{+5.6}$	${9.0}_{-1.2}^{+2.7}$	${8.5}_{-1.5}^{+2.9}$	${10.4}_{-2.5}^{+3.8}$	${10.9}_{-2.0}^{+1.5}$	${9.9}_{-0.7}^{+0.5}$

Note. Each column shows the median and 16th and 84th percentile range over t = 250–450 Myr and t = 250–350 Myr for R8 and LGR4, respectively. Pressure/weight is in units of k_B K cm⁻³, with a multiplication factor of 10³ and 10⁴ for R8 and LGR4, respectively.

Download table as:  ASCIITypeset image

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

4.2.1. 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 $\psi \equiv {G}_{0}\sqrt{T}/{n}_{{\rm{e}}}$ (e.g., Bakes & Tielens 1994; Weingartner & Draine 2001b). The source of FUV radiation is massive young stars (the luminosity-weighted mean age of FUV emitters is ∼10 Myr) so that ${J}_{\mathrm{FUV}}={f}_{\tau }{\dot{S}}_{\mathrm{FUV}}/(4\pi )$ for f_τ , a factor accounting for UV radiative transfer in the dusty ISM, where ${\dot{S}}_{\mathrm{FUV}}={L}_{\mathrm{FUV}}/{L}_{x}{L}_{y}\propto {{\rm{\Sigma }}}_{\mathrm{SFR}}$. ⁹

A simple radiation-transfer solution for uniformly distributed sources in a uniform slab gives

$$\begin{eqnarray}&&{f}_{\tau }=\displaystyle \frac{1-{E}_{2}({\tau }_{\perp }/2)}{{\tau }_{\perp }}\end{eqnarray} \tag{ 22 }$$

(Ostriker et al. 2010). Here, E₂ is the second exponential integral and τ_⊥ = κ_FUVΣ_gas is the mean optical depth to FUV. For κ_FUV = 10³ cm² g⁻¹, f_τ ≈ 1/τ_⊥ at Σ_gas > 20 M_⊙ pc⁻². In the TIGRESS-classic suite, we adopted the approximate form of f_τ as presented in Equation (22) to convert ${\dot{S}}_{\mathrm{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 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)

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({P}_{\mathrm{th},2{\rm{p}}}/{k}_{B})=0.603\,{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{\mathrm{SFR}})+4.99,\end{eqnarray} \tag{ 23a }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{{\rm{\Upsilon }}}_{\mathrm{th}}=-0.506\,{\mathrm{log}}_{10}\left({P}_{\mathrm{DE}}/{k}_{B}\right)+4.45.\end{eqnarray} \tag{ 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 versus Σ_SFR, showing the similar sublinear relationship calibrated to TIGRESS-classic (Equation 23(a); solid black line). Figure 11(d) shows the thermal feedback yield that decreases as P_DE increases (the TIGRESS-classic fit Equation 23(b) 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 a 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. 2023, in preparation).

Zoom In Zoom Out Reset image size

4.2.2. 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⁶ K (or v_SNR ∼ 200 km s⁻¹). The radial momentum injected per SN (p_*) is much larger (∼10⁵ M_⊙ km s⁻¹) than the momentum of the initial SN ejecta (∼10⁴ M_⊙ km s⁻¹) 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 (Kim et al. 2018) and stellar-wind-driven bubbles (Lancaster et al. 2021a, 2021b).

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 a 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 $\rho {v}_{z}^{2}$ 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. 2017a, 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)

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({P}_{\mathrm{turb},2{\rm{p}}}/{k}_{B})=0.96\,{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{\mathrm{SFR}})+6.17,\end{eqnarray} \tag{ 24a }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{{\rm{\Upsilon }}}_{\mathrm{turb}}=-0.06\,{\mathrm{log}}_{10}({P}_{\mathrm{DE}}/{k}_{B})+2.81.\end{eqnarray} \tag{ 24b }$$

Figure 11(b) plots P_turb versus Σ_SFR, showing the expected near-linear relationship (Equation 24(a)). The turbulent feedback yield shown in Figure 11(e) is consistent with the shallow dependence on P_DE seen in Equation 24(b). 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. 2018, 2021).

4.2.3. Magnetic Stress

Given the validity of vertical dynamical equilibrium and agreement of ${ \mathcal 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⁻¹ for model R8-4pc and 720 km s⁻¹ 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 estimated weight:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{\mathrm{SFR}})=1.18\,{\mathrm{log}}_{10}({P}_{\mathrm{tot},2{\rm{p}}}/{k}_{B})-7.43,\end{eqnarray} \tag{ 25a }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{\mathrm{SFR}})=1.17\,{\mathrm{log}}_{10}({ \mathcal W }/{k}_{B})-7.32,\end{eqnarray} \tag{ 25b }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{\mathrm{SFR}})=1.21\,{\mathrm{log}}_{10}({P}_{\mathrm{DE},2{\rm{p}}}/{k}_{B})-7.66.\end{eqnarray} \tag{ 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.

Zoom In Zoom Out Reset image size

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 nonequilibrium 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 TIGRESS-classic 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 Tables 2 and 6, and Figure 13.

Zoom In Zoom Out Reset image size

We summarize the ISM phase distributions by mass and volume within ∣z∣ < 300 pc in Figures 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) occupy 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.

Figures 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 environments. 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. 2023, in preparation), 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⁻¹, and for LGR4, ϒ_tot = 720 km s⁻¹. Both models have similar σ_{eff,2 p} ≈ 12–13 km s⁻¹ and σ_{z,turb,2 p} ∼ 7–8 km s⁻¹.

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 expectations (Ostriker & Kim 2022). The WNM mass fraction may be written as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{M,{\mathtt{WNM}}} & = & \displaystyle \frac{{\rho }_{w}}{{\rho }_{\mathrm{tot}}}{f}_{V,{\mathtt{WNM}}}\\ & = & \displaystyle \frac{{P}_{\mathrm{th}}}{{P}_{\mathrm{tot}}}\displaystyle \frac{{\sigma }_{\mathrm{eff}}^{2}}{{c}_{w}^{2}}{f}_{V,{\mathtt{WNM}}}\end{array}\end{eqnarray} \tag{ 26 }$$

for f_V,WNM the volume-filling factor, where we have used ${P}_{\mathrm{th}}={\rho }_{w}{c}_{w}^{2}$ 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}_{\mathrm{tot}}\equiv {\rho }_{\mathrm{tot}}{\sigma }_{\mathrm{eff}}^{2}$ for ρ_tot, the total midplane density. From Tables 2 and 4, σ_eff and f_V,WNM are very similar between model R8 and LGR4, whereas Table 6 gives a ∼30% lower ratio of thermal-to-total pressure for LGR4 than that for R8. About a 30% decrease in the WNM mass fraction (Table 4 and Figure 13) is consistent with expectation.

5.2.1. Comparison to Milky Way Empirical Constraints on Interstellar Medium Phases

5.2.2. Comparisons to Self-consistent Numerical Models of the Star-forming Interstellar Medium

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₂ conversion factor have been studied in previous work, with Gong et al. (2018, 2020) post-processing the TIGRESS-classic suite with six-ray radiation transfer and steady-state chemistry, and Hu et al. (2021, 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₂ chemistry will be an important tool in understanding observed chemical abundances (see Godard et al. 2022 for CH⁺, and hence warm diffuse H₂, 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 coregulation to the low-metallicity regime, where the thermal feedback yield (and therefore thermal pressure) is expected 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 Section 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) Alfv́en speed, and diffusion by scattering off of MHD waves that are likely self-generated for GeV and lower energies (Armillotta et al. 2021, 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. 2020, 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.

For the photoelectric and Lyman–Werner bands, we calculate the area-averaged intensity at z = z_p-p as a function of the cosine angle $\cos \theta =\hat{{\boldsymbol{k}}}\cdot \hat{{\boldsymbol{z}}}$ (assuming azimuthal symmetry) as

$$\begin{eqnarray}&&\langle I\rangle (\cos {\theta }_{i})\equiv \displaystyle \frac{\int I(x,y,{z}_{{\rm{p}}{\rm{-}}{\rm{p}}};\cos {\theta }_{i}){dxdy}}{\int {dxdy}},\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{1}{{L}_{x}{L}_{y}}\displaystyle \sum _{\mathrm{rays}}\displaystyle \frac{{\rm{\Delta }}{L}_{\mathrm{ray}}}{2\pi \alpha {\rm{\Delta }}x{\rm{\Delta }}\cos \theta },\end{eqnarray} \tag{ A2 }$$

where α = ρ κ_d is the dust absorption cross section per unit volume, Δx is the cell size, ΔL_ray is the amount by which the luminosity of the incoming photon packet is reduced as it passes through a cell at z = z_p-p (Kim et al. 2017b), and the summation is taken over all rays passing through the layer z = z_p-p with $\cos {\theta }_{i}\leqslant \cos \theta \lt \cos {\theta }_{i+1}$. We discretize the cosine angle as $\cos {\theta }_{i}=i{\rm{\Delta }}\cos \theta$ with ${\rm{\Delta }}\cos \theta =2/{N}_{\cos \theta }$ and $i=0,1,...,{N}_{\cos \theta }-1$ and adopt ${N}_{\cos \theta }=64$. Assuming plane-parallel geometry, the radiation energy density and flux in the vertical direction at z > z_p-p are given by

$$\begin{eqnarray}&&\langle {{ \mathcal E }}_{\mathrm{rad},{\rm{p}}{\rm{-}}{\rm{p}}}\rangle =(2\pi {\rm{\Delta }}\cos \theta /c)\displaystyle \sum _{i}\langle I\rangle (\cos {\theta }_{i}){e}^{-{\rm{\Delta }}\tau /\cos {\theta }_{i}},\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}&&\langle {F}_{z,\mathrm{rad},{\rm{p}}{\rm{-}}{\rm{p}}}\rangle =2\pi {\rm{\Delta }}\cos \theta \displaystyle \sum _{i}\langle I\rangle (\cos {\theta }_{i}){e}^{-{\rm{\Delta }}\tau /\cos {\theta }_{i}}\cos {\theta }_{i}\end{eqnarray} \tag{ A4 }$$

where ${\rm{\Delta }}\tau ={\int }_{{z}_{{\rm{p}}{\rm{-}}{\rm{p}}}}^{z}\langle \alpha \rangle {dz}$ is the (area-averaged) dust optical depth integrated from z_p-p to z. Similar calculations are done for z < − z_p-p.

We run two suites of simulations for different convergence tests. In the first suite, we take 10 snapshots (from 200 to 300 Myr) of R8-8pc and LGR4-4pc and post-process the radiation fields, fixing MHD quantities and varying the radiation termination parameters (${d}_{\mathrm{xy},\max }$ and ε_PP). Note that the FUV radiation field at ∣z∣ > z_p-p = 300 pc is set by plane-parallel approximation (see Appendix A.1) rather than the direct radiation transfer. Figure 14 shows the mean radiation energy densities of EUV and FUV at z = 0 and z = L_z /4 with respect to the most accurate results, $({d}_{\mathrm{xy},\max },{\varepsilon }_{\mathrm{PP}})=(4096,0)$ for R8 and (2048, 0) for LGR4. Both EUV and FUV near the midplane are well converged for our fiducial choice (shown as a black star) and generally better converged than those at high-z. The FUV radiation field can be significantly underestimated if ε_PP is large. This implies that low-luminosity sources (photon packets) make a significant contribution to the total FUV radiation field.

Zoom In Zoom Out Reset image size

In the second suite, we restart simulations from Model R8-8pc at 200 Myr with six choices of $({d}_{\mathrm{xy},\max },{\varepsilon }_{\mathrm{PP}})$: Model A (512, 0), Model B (1024, 10⁻⁸), Model C (1024, 0), Model D (2048, 10⁻⁸), Model E (2048, 0), and Model F (4096, 0). Note that Model E is our fiducial choice (identical to R8-8pc). Figure 15 shows the convergence of ISM mass and volume fractions within ∣z∣ < 300 pc, midplane pressures, and feedback yields from this test suite. We present the range of values from the long-term evolution over t ∈ (250, 450) Myr. These quantities are well converged even in Model A with ${d}_{\mathrm{xy},\max }=512\,\mathrm{pc}$ and ε_PP = 0 (blue, leftmost box). Note that we did not repeat the same tests for model LGR4 because the additional expense was not merited, given the good convergence of the midplane radiation field shown in Figure 14.

Zoom In Zoom Out Reset image size

We conclude that our fiducial choices result in good convergence in both radiation field properties and the thermal/dynamical properties of the ISM. It is noteworthy that the accuracy of the radiation transfer can be reduced without affecting too much the thermal and dynamical properties of the ISM (especially, those near the midplane). Unless the radiation field itself is the main subject of the study, our approach with early radiation termination can be used for a large parameter space exploration at a reduced cost.