The GHOSDT Simulations: II. Missing H$_2$ in Simulations of a Self-Regulated Interstellar Medium

Observations in the Galaxy and nearby spirals have established that the HI-to-H$2$ transition at solar metallicity occurs at gas weight of $P{\rm DE}/k_B\approx 10^4 \ \rm K \ cm ^{-3}$, similar to solar neighbourhood conditions. Even so, state-of-the-art models of a self-regulated interstellar medium underproduce the molecular fraction ($R_{\rm mol} \equiv M_{{\rm H}2}/M{HI}$) at solar neighbourhood conditions by a factor of $\approx2-4$. We use the GHOSDT suite of simulations at a mass resolution range of $100-0.25\ M_{\odot}$ (effective spatial resolution range of $\sim 20-0.05\ \rm pc$) run for 500 Myr to show how this problem is affected by modeling choices such as the inclusion of photoionizing radiation, assumed supernova energy, numerical resolution, inclusion of magnetic fields, and including a model for sub-grid clumping. We find that $R_{\rm mol}$ is not converged even at a resolution of 1 $M_{\odot}$, with $R_{\rm mol}$ increasing by a factor of 2 when resolution is improved from 10 to $1\ M_{\odot}$. Models excluding either photoionization or magnetic fields result in a factor 2 reduction in $R_{\rm mol}$. The only model that agrees with the observed value of $R_{\rm mol}$ includes our sub-grid clumping model, which enhances $R_{\rm mol}$ by a factor of $\sim3$ compared with our fiducial model. This increases the time-averaged $R_{\rm mol}$ to $0.25$, in agreement with the Solar circle value, and closer to the observed median value of $0.42$ in regions comparable to the solar neighbourhood in nearby spirals. Our findings show that small-scale clumping in the ISM plays a significant role in H$_2$ formation even in high-resolution numerical simulations.

💡 Research Summary

The paper addresses a long‑standing discrepancy between observations of the atomic‑to‑molecular hydrogen transition in the Milky Way and nearby spirals and the predictions of state‑of‑the‑art self‑regulated interstellar medium (ISM) simulations. Observationally, the HI‑to‑H₂ transition at solar metallicity occurs near a dynamical pressure P_DE/k_B ≈ 10⁴ K cm⁻³, and the molecular‑to‑atomic mass ratio (R_mol = M_H₂/M_HI) in the solar neighbourhood is measured to be ≈0.2–0.3, with a median of ≈0.4 in PHANGS‑ALMA regions that share similar gas and stellar surface densities. However, recent galactic‑patch simulations (e.g., SILC, TIGRESS, zoom‑in studies) typically yield R_mol values of only 0.05–0.2, under‑producing H₂ by a factor of 2–4.

To investigate the origin of this shortfall, the authors employ the GHOSDT (Galaxy Hydrodynamical Simulations of a Supernova‑Driven Turbulent ISM) suite, which uses the GIZMO code with meshless finite‑mass (MFM) hydrodynamics, self‑gravity, and magnetohydrodynamics. The simulations span a mass resolution range from 100 M⊙ down to 0.25 M⊙ (effective spatial resolution 20 pc to 0.05 pc) and are run for 500 Myr in a 1 kpc × 1 kpc × 5 kpc box with Σ_gas = 10 M⊙ pc⁻² and solar metallicity. Star formation is treated on a star‑by‑star basis, with a stochastic IMF sampling that determines each star’s ionizing photon budget and eventual supernova (SN) explosion (standard energy 10⁵¹ erg). The models include a uniform background far‑UV field that scales with the instantaneous star formation rate, on‑the‑fly dust and gas shielding, and a simple treatment of photoionization.

The authors explore five key modeling variations: (1) the fiducial model (including photoionization, magnetic fields, and standard SN energy); (2) omission of photoionization; (3) omission of magnetic fields; (4) varying SN energy by factors of 0.5 and 2; and (5) inclusion of a sub‑grid clumping model. The clumping model assumes a log‑normal density probability distribution below the resolved scale and boosts the effective H₂ formation rate to account for unresolved dense cores. Each variant is simulated for the full 500 Myr, allowing the system to forget its artificial initial conditions and reach a statistically steady state.

Key results: – Resolution dependence: improving mass resolution from 10 M⊙ to 1 M⊙ roughly doubles the time‑averaged R_mol (0.05 → 0.10). Even at 0.25 M⊙ the value only reaches ≈0.12, indicating lack of convergence. – Photoionization and magnetic fields each suppress R_mol by about a factor of two when omitted, confirming that ionizing radiation destroys H₂ and magnetic pressure hampers gas compression. – SN energy variations affect R_mol modestly: doubling the SN energy reduces R_mol by ~30 %, while halving it raises R_mol by ~20 %. – The sub‑grid clumping model has the strongest impact: it raises R_mol by a factor of ≈3, yielding a time‑averaged R_mol ≈0.25, which matches the solar‑circle measurement and lies within the PHANGS‑ALMA median range.

The authors conclude that current galactic‑patch simulations, even at sub‑parsec resolution, miss a substantial fraction of the small‑scale density structure that governs H₂ formation. Including a statistical correction for unresolved clumping is essential to reproduce observed molecular fractions. Moreover, realistic treatment of photoionization and magnetic fields is required; neglecting either leads to significant underestimates of R_mol. The study suggests that future work should combine high‑resolution “zoom‑in” runs with physically motivated sub‑grid models and simultaneously track CO emission and the CO‑to‑H₂ conversion factor to close the gap between theory and observation.