Project description

The THESAN project is a suite of large volume radiation hydrodynamic simulations that self-consistently model the reionization process and the galaxies responsible for it with unprecedented physical fidelity. In fact, from a numerical perspective THESAN provides higher resolution than previous simulations of comparable volume and employs galaxy formation models known to produce physical properties in concordance with observations down to the present-day Universe. Such an approach is ambitious but essential to push the frontier in our understanding of the intergalactic medium and its connection to galaxies during the first billion years after the Big Bang. The following sections include non-specialist and technical introductions to the simulation suite focusing on the unique aspects and central scientific themes of the THESAN project.

Numerical simulations have risen to be a third pillar of research in the last few decades, flanking theoretical and experimental (or observational) efforts. They are particularly relevant in astrophysics because nature prevents scientists from performing controlled experiments due to the large distances, physical sizes, long timescales, or extreme densities and energies involved. In fact, in cosmology the experimental setup cannot be altered as it coincides with part or all of the Universe, a disadvantage that is partially alleviated by the wealth of passive observational data collected by telescopes and detectors across the electromagnetic spectrum. Thus, numerical simulations provide synthetic and controllable versions of the Universe that effectively allow astrophysicists to better understand reality by transcending it. The predictions derived from these virtual experiments are only as good as the physics incorporated into the simulations, which is where the strength and uniqueness of THESAN lies. Specifically, these simulations present novel opportunities to study the early universe by simultaneously including an advanced galaxy formation model, a self-consistent treatment of cosmic reionization, and a model for cosmic dust. Each one of these key ingredients is described in greater detail below.

In the aftermath of the Big Bang, the Universe was comprised of a hot dense plasma of matter and radiation. Cosmic expansion led to rapid gas cooling, with the cosmic microwave background (CMB) forever preserving the moments when electrons and protons recombined to form neutral hydrogen atoms. Throughout its infancy, the dim Universe then went through a process of structure formation induced by the growth of overdensities of dark matter, a still-unknown substance making up approximately 85% of the material in the Universe and only interacting through gravity. Galaxies started to form once rarefied gas could collapse into the central regions of dark matter halos. They then continued to grow and evolve according to the interplay between gravity, hydrodynamical forces, and feedback from stars.

Eventually, within the first billion years, a complex landscape emerged in which various astrophysical phenomena govern the motion and properties of gas, connecting a hierarchy of large to small scale physics throughout the Universe. In particular, countless stars act as cosmic furnaces to not only shine as bright lights that illuminate their surroundings, but also inject large amounts of mass and energy into the interstellar medium (ISM) of galaxies through radiation pressure, stellar winds, and supernovae explosions. Furthermore, across their lifetimes they also fuse lighter elements into heavier ones, which can then form complex molecules like cosmic dust, and even planets. Some of these stars became black holes or other exotic objects at the end of their lives. Finally, nearly every galaxy in the present-day Universe also contains a supermassive black hole that at times may have been active enough to influence the entire galaxy through their powerful energy production.

All of these processes are carefully included in the THESAN simulations. The gravitational and hydrodynamical forces are accurately calculated using the novel moving mesh approach with the AREPO code. The entire simulation volume is divided into approximately ten billion individual volume elements, which roughly follow the motion of the gas. In this way, they naturally cluster in regions of high density, where the resolution is needed the most as these are the epicenters of galaxy formation. The other processes mentioned above are implemented following the largely successful IllustrisTNG galaxy formation model, that was developed to produce realistic galaxy properties and populations in broad agreement with available observations of the local and high-redshift Universe. To make all of this possible, the simulations were run on one of the largest supercomputers in the world, the SuperMUC-NG machine at the LRZ (in Germany), employing nearly 30 million hours of continuous calculations distributed across almost 60,000 processors working together simultaneously.

Once the first stars — and subsequently galaxies — formed in the infant Universe, they began to pour large amounts of high-energy photons into the surrounding intergalactic medium. This gradually transformed the ambient neutral hydrogen gas into a hot highly-ionized plasma, a process known as cosmic reionization. Reionization represents the last phase transformation that the Universe went through and forever changed the stage upon which galaxies formed and evolved. Unfortunately, it is extremely challenging to gather data on this epoch as the large amount of time between us and when this transformation took place translates to truly astronomical distances. However, a detailed chronology and characterization of cosmic reionization will continue to emerge in the coming years.

The THESAN simulations are specifically designed to advance our knowledge regarding cosmic reionization. For this reason, they combine the realistic IllustrisTNG galaxy formation model with an accurate radiation transport algorithm to track the propagation of light from stars throughout the Universe. Given the colossal speed of light compared to typical gas motions, modelling radiation in dynamical simulations is one of the most challenging problems in computational astrophysics. However, this approach unlocks a new level of physical fidelity that is unsurpassed by other methods due to the combination of scales, resolution, and richness of physical processes explored. The THESAN simulations thus enable us to simultaneously investigate the formation of the first galaxies and their impact on the diffuse gas filling the space between them. Thus, THESAN provides an exquisite and invaluable tool for state-of-the-art scientific research, facilitating the interpretation of forthcoming data from current and upcoming telescopes.

The final novel ingredient of the THESAN simulations is the inclusion of cosmic dust, which renders this project even more unique within the landscape of astrophysical simulations. Cosmic dust indicates an ensemble of tiny grains of solid material that are formed in galaxies and are found between stars. One side effect is that dust locks a fraction of the atomic and molecular gases made from heavy elements into amorphous and crystalline structures, which can have an important impact on cosmic chemistry and cooling. Dust also transforms ultraviolet and optical light into lower-energy photons, thereby changing how radiation influences galaxy formation processes and in some cases allowing galaxies to shine brightly at infrared wavelengths. Thus, dust reprocessed light has recently become a powerful method for probing the properties of the first galaxies, and is rapidly providing a wealth of information about the conditions of the infant Universe. For this reason, although dust modelling in simulations is relatively new territory, self-consistently simulating this component is of primary importance to test our knowledge of the first galaxies with upcoming observations from increasingly powerful telescopes.

Gaining insights into the Epoch of Reionization has proven to be quite challenging because even the deepest Hubble Space Telescope (HST) observations are only able to detect galaxies during the late stages of reionization. Observations of high redshift quasars have provided insights into the state of the large scale intergalactic medium (IGM), as well as important constraints on the duration of reionization. Unfortunately, it is quite difficult to build an overall picture of the state of the IGM using the information coming from individual sight-lines to quasars alone. However, the state of the field is about to change due to a suite of instruments and telescopes that will become available in the next few years. For example, the James Webb Space Telescope (JWST) with its large mirror and infrared frequency coverage, will enable the discovery of fainter and higher-redshift galaxies and constrain their rest-UV/optical emission. Additionally, the fragmented view of the IGM will soon be replaced by a more complete picture of the neutral hydrogen distribution in the Universe by current and upcoming 21 cm line intensity mapping (LIM) experiments like the Low-Frequency Array (LOFAR), Hydrogen Epoch of Reionization Array (HERA), and Square Kilometre Array (SKA). In a similar vein, there are plans to use LIM of spectral lines originating from many individually unresolved galaxies to faithfully trace the growth and evolution of large-scale structures. These advancements will unleash a flood of high-redshift observational data that promises to usher in a new era of cosmic reionization studies. The THESAN simulations have been specifically designed to fully exploit the upcoming observations by making accurate predictions of the large-scale reionization process and resolved properties of the sources responsible for it.

We outline below a few of the key science areas that motivated the THESAN simulations:

• High redshift galaxy formation and JWST predictions: JWST will be launched this year and will transform our understanding of the formation and evolution of the first galaxies. Presently, the galaxy population at z>4 is essentially unexplored, with HST discovering just about thousand galaxy candidates at z=6-8 and only a few at higher redshifts. Most of these galaxies do not have a spectroscopic confirmation/characterization and only a handful of objects have been observed with ALMA. Hence, current observational constraints only provide weak tests of reionization-era galaxy-formation models. JWST will change the game thanks to its imaging and spectroscopic capabilities in the infrared that will enable the discovery of fainter, higher-redshift galaxies (up to z~20) and constrain their rest-UV/optical emission. This will provide measurements of, among other things, (i) spectroscopic redshifts; (ii) accurate stellar masses by breaking degeneracies between age, dust and nebular emission; (iii) galaxy structure, including the formation of the first disks and the appearance of the first bulges; (iv) emission-line properties such as ionization (AGN/shocks), metallicity and AGN/star-formation driven outflows. THESAN provides crucial, realistic predictions for all these observations. The incorporation of a successful galaxy formation model allows robust forecasts for future JWST surveys. These, in turn, will enable us to test our knowledge of processes involved in galaxy formation out to an unprecedented redshift. This is of paramount importance in the understanding of how galaxies emerged from the initial, tiny fluctuations of the density field, as well as of how they affected their environment.
• Using the Lyα forest to study the IGM: By taking advantage of bright quasar sources as backlights, the physical state of the IGM can be studied through the Lyman-alpha (Lyα) forest. Historically, this technique has been limited to the post-reionization Universe as a consequence of the very-large cross-section of the Lyα transition. However, in recent years improvements in instrumentation have allowed us to push this kind of observations deep into the tail-end of reionization, unveiling features that are largely unexplored. This is a promising ground for advancing our knowledge of the Epoch of Reionization (EoR). THESAN allows us to investigate the inception of transmission regions in the Lyα forest. These features are produced by single, bright galaxies and therefore require the inclusion of a state-of-the-art galaxy formation model as we use here. The rarity of these features at high-redshift and the patchy nature of the reionization process demands large, simulated volumes, in order to properly capture statistical variations in the IGM properties. Additionally, recent observations have shown that the IGM optical depth distribution is difficult to reproduce in the standard EoR model. This is especially interesting as this is a first step beyond globally averaged measures of the IGM properties. Many different mechanisms have been suggested to reconcile the observed and predicted distributions, but no conclusive answer has been agreed upon by the scientific community. THESAN allows us to investigate this further and provide a robust prediction from a suite of full-physics state-of-the-art simulations.
• Predictions for line intensity mapping experiments: Observations of the hyperfine transition of neutral hydrogen (best known as the 21 cm line) promises to revolutionize our understanding of the early Universe. It will allow us to transition from the discrete view provided by single lines of sight towards bright objects to a tomographic picture. This is the target of current and upcoming instruments like the SKA and HERA. The fluctuations in the 21 cm map will depend on a range of physical processes inside and around galaxies, such as feedback, recombinations, escape fractions, etc. Additionally, LIM experiments such as COMAP, CCAT-p and SPHEREx aim to map the molecular and nebular emission lines emerging from the ISM of these early galaxies. The intensity and prevalence of these lines will depend substantially on the properties of the gas in high redshift galaxies. LIM is advantageous to galaxy surveys as it is sensitive to all sources of emission in the line and thus enables the universal study of galaxy formation and evolution. Intensity maps of [CII], [OI], and [OIII] fine structure lines, CO rotational transitions, and Hα emission, among others, may be used to trace the galaxy distribution in the same cosmological volume as the 21 cm observations. The novelty of the THESAN simulations, which combines a state-of-the-art galaxy formation model with RHD and at the same time simulate large volumes allows for predictions of both the nebular lines and 21 cm emission in a consistent manner. Using both these observations individually and through cross-correlations will provide important insights into both the sources of reionization and how they combine with the IGM to ionize it.
• Lyα emission and transmission: The Lyα line of atomic hydrogen is among the strongest lines emitted by high-redshift galaxies, which makes it an excellent target for spectroscopic redshift determinations. Observational campaigns with the JWST will provide unprecedented samples of Lyα emitting galaxies (LAEs) that will serve as probes of reionization and reveal properties about star-forming galaxies. The THESAN and forthcoming THESAN-ZOOMS projects allow us to self-consistently connect the intrinsic Lyα emission from recombinations and collisional excitation, resonant scattering within the interstellar and circumgalactic medium (ISM/CGM), and frequency-dependent transmission through the IGM. Accurate Monte Carlo radiative transfer calculations for ionizing, Lyα, Hα, continuum, and other forms of radiation reveal the physical mechanisms at play in high-z LAEs and test if their properties are similar to those identified at lower redshifts, e.g. in terms of dust content and star formation. THESAN is ideal for exploring why galaxies become LAEs, signatures for Lyα line intensity mapping, factors regulating IGM transmissivity, line-of-sight escape fractions, Lyα halos and their connection to the CGM, and links between spatially/spectrally resolved properties and galaxy formation and evolution. By comparing Lyα simulations and observations we can extract additional missing information about high-z galaxies to help guide future studies of the EoR.

• Radiation transport: The propagation of radiation is handled using a moment-based approach to solve the radiation transport equations. We solve the set of coupled hyperbolic conservation equations for photon number density and photon flux. This set of equations is closed using the M1 scheme. The Riemann problem is solved at each cell interface by computing the flux using Godunov’s approach. We achieve second order accuracy by replacing the piecewise constant (PC) approximation of Godunov’s scheme with a slope-limited piece-wise linear spatial extrapolation and a half timestep, first order time extrapolation, obtaining the primitive variables on both sides of the interface. This approach is fully local and does not scale with the number of sources in the simulation box, a feature particularly appealing given the large amount of stellar particles present in the simulations.
• IllustrisTNG galaxy formation model: The simulations employ the state-of-the-art IllustrisTNG galaxy formation model, which updates the previous Illustris simulations with a revised kinetic AGN feedback model for the low accretion state and an improved parameterization of galactic winds. It includes: (i) a sub-resolution treatment of the interstellar medium (ISM) as a two-phase gas where cold clumps are embedded in a smooth, hot phase produced by supernova explosions; (ii) feedback from supernova explosions and stellar winds, in the form of kinetic and thermal energy; (iii) the production and evolution of nine elements (H, He, C, N, O, Ne, Mg, Si and Fe), as well as the tracking of the overall gas metallicity; (iv) density-, redshift-, metallicity- and temperature-dependent cooling; and (v) black-hole formation (via a seeding prescription), growth and feed-back in two different regimes (quasar- and radio-mode). All of these elements are key to realistically simulating the reionization history of the Universe, which heavily depends on the source properties, and to investigating the correlation between small and large scales in the Universe.
• An empirical dust model: The IllustrisTNG galaxy formation model is augmented with a set of empirical relations that describe the production and destruction of cosmic dust. We follow the approach of McKinnon et al. (2016, 2017) and numerically treat dust as a property of the gas resolution elements. Within individual resolution elements the model tracks the mass of dust in five chemical species (C, O, Mg, Si and Fe). The dust mass can increase via two processes: dust production during stellar evolution, and dust growth via deposition of gas-phase metals and decrease due to sputtering in high temperature gas.
• Simulations: THESAN is a suite of large volume(L_box=95.5 cMpc) radiation hydrodynamic simulations that self-consistently model the hydrogen reionization process and the sources responsible for it. The highest resolution simulation THESAN-1 has a DM mass resolution of 3.12 x 10⁶M_☉ and a baryonic mass resolution of 5.82 x 10⁶M_☉. The gravitational forces are softened on scales of 2.2 ckpc with the smallest cell radii reaching 10 pc. This allows us to model atomic cooling halos throughout the entire simulation volume. A suite of lower resolution (8 times lower mass and 2 times lower spatial resolution) simulations helps quantify the changes to reionization induced by halo mass dependent escape fractions (THESAN-HIGH-2 and THESAN-LOW-2), alternative dark matter models (THESAN-SDAO-2), assuming an instantaneous reionization model (THESAN-TNG-2 and THESAN-TNG-SDAO-2), back reaction of the reionization process on galaxy formation (THESAN-NORT-2), and numerical convergence (THESAN-2 and THESAN-WC-2). The full THESAN simulation set and the associated parameters are cataloged in the following table:

Figure: Global properties of the IGM in the THESAN simulations, compared to observations. From left to right: volume-averaged hydrogen neutral fraction evolution, volume-averaged temperature at mean density evolution, and CMB Thomson optical depth as a function of maximum integration redshift. The bottom panels show the variance around the THESAN-1 curve, computed by subdividing the original box into 8 (inner shaded region) and 64 (outer shaded region) sub-boxes. The THESAN simulations agree well with available data.

We highlight here a few of the early science results from the first THESAN papers:

• Reionization histories: The fiducial simulation and model variations all produce realistic reionization histories and match the observed neutral fraction evolution and the optical depth to the CMB. The fiducial simulation has relatively large neutral regions below z<6, mimicking ‘late’ reionization models, which have recently been invoked to explain long troughs in the Lyα forest.
• High-redshift galaxies: The simulations produce a realistic high redshift galaxy population. They match a wealth of high-redshift observationally inferred data, including the stellar-to-halo-mass relation, galaxy stellar mass function, dust-attenuated ultra-violet luminosity functions, and the mass-metallicity relation, despite the galaxy formation model being mainly calibrated at z=0.
• IGM-galaxy connections: The modulation of the transmitted Lyα flux as a function of the distance from galaxies is a particularly stringent test, since it depends simultaneously on the properties of reionization and of the galaxy population. We find that the simulations can match the general observed trend, only if the reionization completes below z<5.2.
• Lyα IGM transmission: The simulations show that large ionized bubbles form around the brightest galaxies at any epoch. Such accelerated reionization regions provide an advantage for Lyα transmission of red peaks, thereby boosting the visibility of LAEs around UV bright galaxies. However, we also show that massive galaxies have larger covering fractions near line center as infall motion and self-shielding structures play significant roles in shaping the IGM transmissivity.
• Predictions for 21 cm measurements: The simulations predict that different reionization models produce different bubble size distributions at the same ionization fraction. For example, reionization scenarios where low mass halos dominate the ionizing output form a large number of very small bubbles, while models in which high mass halos contribute the most number ionizing photons form less numerous but large HII regions. Encouragingly, these different bubble sizes show unique signatures in the 21 cm power spectrum especially on the slope of the power spectrum at large spatial scales, enabling current and upcoming 21 cm experiments to accurately characterize the sources that dominate the ionizing photon budget.

Visit the publications page for more science results utilizing the THESAN simulations.