THESAN-1: Flagship Reionization Simulation

• THESAN-1 is a flagship simulation modeling the Epoch of Reionization by coupling the IllustrisTNG galaxy-formation model with on-the-fly radiative transfer in a 95.5 cMpc volume.
• Its high resolution resolves atomic-cooling haloes and accurately captures the evolution of ionized bubbles, galaxy clustering, and ionizing escape fractions critical for reionization studies.
• The simulation underpins extensive analyses—including IGM thermal structure, Lyman-α emitters, and bubble merger histories—while serving as a benchmark for derived zoom and dark-matter-only programs.

THESAN-1 is the flagship simulation of the THESAN project, a large-volume radiation-hydrodynamic program designed to model the Epoch of Reionization by evolving galaxy formation and ionizing-radiation transport in a single framework. It couples the IllustrisTNG galaxy-formation model to on-the-fly radiative transfer in a $95.5\,\mathrm{cMpc}$ cosmological volume, with sufficient resolution to resolve atomic-cooling haloes while retaining the large-scale ionized-bubble network required for reionization studies. As a result, THESAN-1 has become a common reference simulation for the high-redshift intergalactic medium, ionized bubble growth, ionizing escape fractions, Lyman-$\alpha$ transmission, Lyman-$\alpha$ emitter statistics, and several derived programs including zoom simulations and Local Group reionization mappings (Kannan et al., 2021, Garaldi et al., 2021).

1. Position within the THESAN program

THESAN is a suite of cosmological radiation-magnetohydrodynamic or radiation-hydrodynamic simulations built to bridge galaxy formation and reionization studies in a fully coupled setting. Within that suite, THESAN-1 is the highest-resolution large-volume full-physics run and serves as the fiducial reference against which medium-resolution source-variation, convergence, and alternative-dark-matter runs are compared. Its scientific role is twofold. First, it provides a representative reionization volume large enough to capture bubble statistics, galaxy clustering, and late-stage percolation. Second, it resolves the galaxy population down to the atomic-cooling threshold, which is required to model the sources that dominate the ionizing photon budget during much of the Epoch of Reionization (Kannan et al., 2021).

The same simulation underlies several later THESAN analyses. It is the basis for studies of IGM thermal and ionization structure, effective Lyman-$\alpha$ optical-depth fluctuations, ionized-bubble environments, escape fractions, Lyman-$\alpha$ luminosity-function calibration, line-of-sight bubble statistics around Lyman-$\alpha$ emitters, and bubble merger histories. In that sense, “THESAN-1” denotes not only a specific numerical realization but also the common physical background on which a large fraction of the THESAN literature is constructed.

2. Numerical realization and physical model

THESAN-1 evolves a periodic box of comoving side length $L_{\rm box}=95.5\,\mathrm{cMpc}$ with $2100^3$ dark matter particles and $2100^3$ initial gas cells. The dark matter particle mass is $m_{\rm DM}=3.12\times 10^{6}\,M_\odot$, the initial gas cell mass is $\alpha$0, the gravitational softening for collisionless components is $\alpha$1, and minimum gas cell sizes reach $\alpha$2 at $\alpha$3. The run is evolved down to $\alpha$4, and this resolution resolves atomic-cooling haloes across the full volume (Kannan et al., 2021, Xu et al., 2022).

Property	THESAN-1 value	Note
Box size	$\alpha$5	Periodic volume
Resolution elements	$\alpha$6 DM + $\alpha$7 gas	Flagship run
Dark matter mass	$\alpha$8	Per particle
Gas mass	$\alpha$9	Initial cell mass
Softening	$\alpha$0	Collisionless components
Minimum cell size	$\alpha$1	At $\alpha$2
End redshift	$\alpha$3	Full EoR coverage

The code base is AREPO-RT, the radiation-hydrodynamic extension of AREPO. Gas dynamics are solved on a moving Voronoi mesh, gravity uses a Tree-PM scheme, and ionizing radiation is evolved with a moment-based transport method using M1 closure. The radiation field is discretized by ionization thresholds at 13.6, 24.6, 54.4, and $\alpha$4 eV, and a reduced speed of light $\alpha$5 is adopted. Non-equilibrium thermochemistry tracks H/He ionization states and heating/cooling. The galaxy-formation model is IllustrisTNG, including star formation with a Springel-Hernquist multiphase ISM, stellar feedback, black-hole growth, AGN feedback, metal enrichment, and a dust model following McKinnon et al. Stellar ionizing spectra are taken from BPASS with a Chabrier IMF, and THESAN-1 uses a sub-resolution birth-cloud escape fraction $\alpha$6 to match reionization constraints (Kannan et al., 2021, Garaldi et al., 2021).

This architecture gives THESAN-1 a distinctive niche. It is large enough to contain rare, bright galaxies and large ionized regions, but fine-grained enough to model the faint-galaxy population and the small-scale absorbers that regulate ionizing-photon transport. That combination explains why later THESAN work repeatedly adopts THESAN-1 as the backbone simulation for post-processing studies.

3. Reionization history, IGM structure, and ionized-bubble topology

THESAN-1 produces a realistic late reionization history that matches the observed evolution of the global neutral hydrogen fraction and the electron-scattering optical depth, while also reproducing the rapid evolution of the ionizing mean free path at high redshift (Kannan et al., 2021, Garaldi et al., 2021). Its synthetic Lyman-$\alpha$7 forest spectra yield effective optical-depth evolution consistent with available data, although the distribution suggests that an even-later reionization than simulated may still be favored. The simulation also reproduces, for the first time in this program, the modulation of Lyman-$\alpha$8 flux as a function of galaxy distance, and later work showed that the effective optical depth of quasar sightlines is most sensitive to galaxies at a redshift-dependent distance that is $\alpha$9 at the tail end of reionization (Garaldi et al., 2024).

A major THESAN-1 development was the construction of local bubble-size fields with the mean-free-path method. In this framework, each point in the ionized IGM is assigned an effective bubble size $\alpha$0, enabling volume-weighted bubble-size distributions, environment-conditioned statistics, and direct comparison to Lyman-$\alpha$1 and 21-cm observables. THESAN-1 shows that regions ionized early undergo a slow growth period, whereas regions ionized later often experience a rapid “flash ionization” when they are engulfed by a large pre-existing bubble. High-overdensity regions have larger characteristic bubble sizes, but this correlation weakens as reionization progresses and large percolated bubbles come to dominate the volume (Neyer et al., 2023).

A complementary analysis based on the spatially resolved redshift of reionization $\alpha$2 and a bubble merger-tree construction identifies three stages of bubble evolution. The first is initial slow expansion around the earliest ionizing sources. The second is accelerated growth through percolation as bubbles begin to merge. The third is rapid expansion dominated by the largest bubble. In THESAN-1, the largest bubble emerges by $\alpha$3–10, well before the midpoint of reionization, and later expands to encompass almost the entire ionized volume. The same study finds a sharp decline in the number of bubbles with radii around $\alpha$4 compared to smaller sizes, indicating a characteristic scale in the final segmented bubble-size distribution (Jamieson et al., 2024).

These results together imply a strongly environment-biased, percolative, inside-out reionization morphology. They also clarify why simple spherical-bubble descriptions become increasingly inaccurate once the bubble network enters the percolation regime.

4. Escape fractions and the source population of reionization

A central THESAN-1 application is the calculation of ionizing escape fractions for reionization-era galaxies. In post-processing, the halo escape fraction is defined as

$\alpha$5

where $\alpha$6 is the intrinsic ionizing photon production rate inside $\alpha$7 and $\alpha$8 is the rate crossing the virial sphere into the IGM (Yeh et al., 2022).

THESAN-1 finds that the escape fraction depends strongly on halo mass and redshift. Low-mass galaxies with $\alpha$9 are the main drivers of reionization above $\alpha$0, while high-mass galaxies with $\alpha$1 dominate the escaped ionizing photon budget at lower redshifts. The halo-to-halo variation in escape fraction decreases for higher-mass haloes, and this is associated with more settled galactic structure, greater SFR stability, and a larger fraction of sightlines within each halo contributing significantly to the escaped flux. Dust is capable of reducing the escape fractions of massive galaxies, though the impact on the global $\alpha$2 depends on the dust model. AGN are unimportant for reionization in THESAN, and their escape fractions are lower than stellar ones because they are located near the centers of galaxy gravitational potential wells (Yeh et al., 2022).

This division of labor between low-mass and high-mass galaxies is consistent with the broader THESAN-1 reionization narrative: faint systems initiate and sustain early ionization growth, whereas more massive galaxies, which are statistically underrepresented in many smaller-box simulations, become essential for the rapid late stages of reionization.

5. Lyman-$\alpha$3, Lyman-$\alpha$4 emitters, and line-of-sight bubble statistics

THESAN-1 has also been used to construct empirical but physically anchored models of Lyman-$\alpha$5 emitter observables. An early THESAN analysis calibrated a five-parameter dust-and-spectral model to observed Lyman-$\alpha$6 luminosity functions at $\alpha$7 and $\alpha$8, using Gaussian Process Regression with SWIFTEmulator. In that framework, intrinsic Lyman-$\alpha$9 luminosities are taken directly from resolved recombinations, collisional excitation, and unresolved H II regions, while IGM transmission is computed by ray tracing through the native THESAN-1 mesh (Xu et al., 2022). The intrinsic production rate is written as

$\alpha$0

A later THESAN-1 study replaced that effective treatment with a six-parameter LAE model combining intrinsic Lyman-$\alpha$1 production, an outflow-shaped emergent spectral profile, dust attenuation tied to UV continuum extinction, and exact line-of-sight IGM transmission from 768 rays per central galaxy. The final observed luminosity is

$\alpha$2

Using this framework, THESAN-1 shows that before the midpoint of reionization, galaxies in larger line-of-sight bubbles, $\alpha$3, have higher observed Lyman-$\alpha$4 luminosity and equivalent width, and that these correlations weaken as percolation progresses and the IGM becomes increasingly ionized. In LAE-selected samples with $\alpha$5, Lyman-$\alpha$6 properties correlate with bubble size more strongly than UV magnitude, especially at $\alpha$7. The evolving LAE fraction is therefore a practical statistical tracer for bubble size (Neyer et al., 21 Oct 2025).

THESAN-1 now supplies a mapping between LAE selections and bubble-size statistics, and public catalogues are planned to include, for each galaxy and each of 768 sightlines, intrinsic and observed Lyman-$\alpha$8 luminosities, $\alpha$9, IGM transmission curves and integrated $L_{\rm box}=95.5\,\mathrm{cMpc}$0, UV magnitudes, equivalent widths, LOS bubble sizes $L_{\rm box}=95.5\,\mathrm{cMpc}$1, and environmental overdensities (Neyer et al., 21 Oct 2025).

6. Derived programs, public legacy, and limitations

THESAN-1 functions as the parent environment for several derived programs. THESAN-ZOOM adopts regions selected from the THESAN volume and imports the large-scale radiation field as a boundary condition, enabling studies of external UV suppression, bursty star formation, and lingering Population III star formation in a realistic patchy reionization context (Zier et al., 4 Mar 2025, McClymont et al., 28 Feb 2025, Zier et al., 5 Mar 2025). THESAN-DARK-1, the dark-matter-only twin, is cross-matched to THESAN-1 to map reionization histories onto $L_{\rm box}=95.5\,\mathrm{cMpc}$2 Local Group analogues, showing that overdensity and proximity to Virgo-like protoclusters strongly modulate $L_{\rm box}=95.5\,\mathrm{cMpc}$3 and can generate pairwise reionization offsets of up to $L_{\rm box}=95.5\,\mathrm{cMpc}$4 (Zhao et al., 22 Jul 2025). THESAN-1 has also been used to make first-generation Lyman-$L_{\rm box}=95.5\,\mathrm{cMpc}$5 intensity-mapping forecasts, where emission-only power lies above SPHEREx sensitivity but absorption-included predictions are $L_{\rm box}=95.5\,\mathrm{cMpc}$6 orders of magnitude lower, emphasizing the importance of resonant scattering and subgrid galactic radiative transfer (Almualla et al., 5 Dec 2025).

Its limitations are equally well documented. The finite volume underrepresents the rarest, brightest galaxies and the most extreme opaque sightlines, and it can underproduce the highest-luminosity bins of the Lyman-$L_{\rm box}=95.5\,\mathrm{cMpc}$7 luminosity function (Neyer et al., 21 Oct 2025, Garaldi et al., 2024). Reionization in THESAN-1 is slightly early relative to some of the latest forest-based constraints (Almualla et al., 5 Dec 2025). ISM-scale Lyman-$L_{\rm box}=95.5\,\mathrm{cMpc}$8 radiative transfer remains unresolved, so emergent line profiles and dust attenuation are treated with calibrated effective models rather than first-principles Monte Carlo transport through a multiphase ISM (Xu et al., 2022, Neyer et al., 21 Oct 2025). For intensity mapping, scattering back into the line of sight is neglected in the most conservative calculations, implying that some THESAN-1 LIM predictions are lower limits rather than full observables (Almualla et al., 5 Dec 2025).

Even with those caveats, THESAN-1 occupies a distinctive place in reionization studies. It is simultaneously a large-scale ionization simulation, a resolved high-redshift galaxy simulation, a platform for Lyman-$L_{\rm box}=95.5\,\mathrm{cMpc}$9 and 21-cm forward modeling, and the parent environment for several zoom and descendant-mapping projects. Its importance lies less in any single metric than in the coherence of that framework: galaxy formation, ionizing escape, bubble growth, IGM transmission, and observable selection effects are all represented within one common numerical realization.