[Paper Review] Introducing the THESAN project: radiation-magnetohydrodynamic simulations of the epoch of reionization

Rahul Kannan, Enrico Garaldi+5 more|arXiv (Cornell University)|Oct 1, 2021

Galaxies: Formation, Evolution, Phenomena6 references13 citations

TL;DR

The THESAN project presents a large-volume (95.5 cMpc), radiation-magnetohydrodynamic simulation suite that self-consistently models galaxy formation and intergalactic medium (IGM) ionization during the Epoch of Reionization (EoR). Using the arepo-rt radiation hydrodynamics solver and the IllustrisTNG galaxy formation model with non-equilibrium thermochemistry and dust, it reproduces observed reionization histories and high-redshift galaxy properties, demonstrating that distinct bubble size distributions from different ionizing source populations imprint unique signatures on the 21 cm power spectrum, especially at large spatial scales.

ABSTRACT

We introduce the THESAN project, a suite of large volume (L = 95.5 cMpc) radiation-magnetohydrodynamic simulations that simultaneously model the large-scale statistical properties of the IGM during reionization and the resolved characteristics of the galaxies responsible for it. The flagship simulation has dark matter (DM) and baryonic mass resolutions of $3.1 imes 10^6$ $M_\odot$ and $5.8 imes 10^5 $ $M_\odot$, respectively. The gravitational forces are softened on scales of 2.2 ckpc with the smallest cell sizes reaching 10 pc at z=5.5, enabling predictions down to the atomic cooling limit. The simulations use an efficient radiation hydrodynamics solver (AREPO-RT) that precisely captures the interaction between ionizing photons and gas, coupled to well-tested galaxy formation (IllustrisTNG) and dust models to accurately predict the properties of galaxies. Through a complementary set of medium resolution simulations we investigate the changes to reionization introduced by different assumptions for ionizing escape fractions, varying DM models, and numerical convergence. The simulations produce realistic reionization histories that match the observed evolution of the global neutral hydrogen fraction and electron scattering optical depth to reionization. They also match a wealth of high-redshift observationally inferred data, including the stellar-to-halo-mass relation, stellar mass function, star formation rate density, and the mass-metallicity relation, despite the galaxy formation model being mainly calibrated at z=0. We demonstrate that different reionization models give rise to varied bubble size distributions that imprint unique signatures on the 21cm emission, especially on the slope of the power spectrum at large spatial scales, enabling current and upcoming 21cm experiments to accurately characterise the sources that dominate the ionizing photon budget. [abridged]

Motivation & Objective

• To model the coupled evolution of galaxies and the ionized intergalactic medium (IGM) during the Epoch of Reionization (EoR) with high resolution and physical fidelity.
• To produce a self-consistent simulation suite that matches global reionization observables, including the neutral hydrogen fraction and electron scattering optical depth.
• To investigate how different assumptions about ionizing escape fractions, dark matter models, and numerical resolution affect reionization history and IGM morphology.
• To enable predictions and interpretation of upcoming observations from JWST, HERA, SKA, and other high-redshift telescopes.
• To demonstrate that 21 cm power spectrum slopes at large scales can distinguish between reionization models dominated by low-mass versus high-mass halos.

Proposed method

• The simulations use the arepo-rt radiation hydrodynamics solver to accurately model ionizing photon propagation and gas heating in the IGM.
• The galaxy formation model is based on IllustrisTNG but enhanced with non-equilibrium thermochemistry and an empirical dust model to improve high-redshift physics.
• The fiducial simulation has a box size of 95.5 cMpc, with dark matter and baryonic mass resolutions of 3.1×10⁶ M☉ and 5.8×10⁵ M☉, respectively.
• Gravitational softening is set to 2.2 ckpc, with spatial resolution reaching 10 pc at z=5.5, enabling resolution down to the atomic cooling limit.
• A suite of medium-resolution simulations explores variations in ionizing escape fractions, dark matter models, and numerical convergence.
• The simulations are calibrated to match observational constraints on the global neutral hydrogen fraction, electron scattering optical depth, stellar-to-halo mass relation, and galaxy stellar mass function.

Experimental results

Research questions

• RQ1 How do different assumptions about ionizing escape fractions affect the reionization history and IGM morphology?
• RQ2 What is the impact of varying dark matter models on the evolution of H II regions and reionization timing?
• RQ3 How do different galaxy formation models, particularly those with varying ionizing source populations, affect the 21 cm power spectrum at large spatial scales?
• RQ4 To what extent can the 21 cm power spectrum slope at small k distinguish between reionization models dominated by low-mass versus high-mass halos?
• RQ5 Can the THESAN simulations reproduce observed high-redshift galaxy properties such as the stellar mass function and mass-metallicity relation despite being calibrated only at z=0?

Key findings

• The fiducial THESAN simulation successfully reproduces the observed evolution of the global neutral hydrogen fraction and electron scattering optical depth to reionization.
• The simulations match high-redshift observational constraints, including the stellar-to-halo mass relation, galaxy stellar mass function, star formation rate density, and mass-metallicity relation, despite galaxy formation model calibration at z=0.
• Different reionization models—driven by low-mass or high-mass halos—produce distinct bubble size distributions in the IGM.
• These differing bubble size distributions leave unique imprints on the 21 cm power spectrum, particularly in the slope at large spatial scales (k ≈ 0.2 h/cMpc).
• The slope of the 21 cm power spectrum at small k is shown to be a robust discriminator between reionization models, enabling future constraints on the dominant ionizing sources.
• The simulations demonstrate that high-resolution radiation-MHD modeling is essential for predicting accurate IGM and galaxy properties during the EoR, with implications for interpreting data from JWST, HERA, and SKA.