Cosmological RHD Simulations

• Cosmological RHD simulations are computational models that couple dark matter, gas dynamics, radiative transfer, and non-equilibrium chemistry to study cosmic reionization and galaxy formation.
• They employ advanced numerical methods, including AMR, SPH, and GPU-accelerated solvers, to resolve both large-scale structures and fine-scale radiation–hydrodynamics interactions.
• These simulations provide critical insights into reionization timing, feedback effects on low-mass galaxies, and observable signatures like the 21 cm line and UV luminosity functions.

Cosmological radiation-hydrodynamic (RHD) simulations are computational models that jointly solve for the evolution of gas, dark matter, radiation fields, and typically non-equilibrium chemistry, across cosmological volumes. They enable direct investigation of the interplay between radiative feedback, star/galaxy formation, intergalactic medium (IGM) reionization, and observable tracers (e.g., 21 cm line, UV luminosity functions, kSZ, metal lines) from the dawn of structure formation through to later cosmic epochs. These simulations are essential for interpreting and predicting the observable signatures of cosmic reionization, the thermal/ionization state of the IGM and circum/interstellar media, and the impact of radiative processes on galaxy formation.

1. Physical and Mathematical Foundations

Cosmological RHD simulations solve a coupled set of equations describing:

• Collisionless dark matter: Typically advanced with a particle-mesh (PM) or tree/particle-particle method, evolving the Vlasov–Poisson system to capture structure formation at large scale.
• Gas dynamics: Euler equations in comoving coordinates for mass, momentum, and energy, sometimes formulated on a moving mesh or using smoothed-particle hydrodynamics (SPH; e.g., GADGET, Gasoline2) or adaptive mesh refinement (AMR; e.g., RAMSES, Enzo).
• Radiative transfer (RT): Typically formulated as either angle-resolved RT:

$$\frac{1}{c}\frac{\partial I_\nu}{\partial t} + \mathbf{n}\cdot\nabla I_\nu = j_\nu - \kappa_\nu I_\nu$$

or via moment equations with a closure (e.g., M1) for computational tractability.

• Thermochemistry: Non-equilibrium rate equations track hydrogen and helium ionization states, with photoionization, recombination, and explicit radiative heating/cooling, sometimes extended to metal line cooling and chemistry (e.g., Aurora, Enzo runs).
• Star formation and feedback: Subgrid models for star formation, supernova (SN) feedback, and often radiative escape fractions ($f_{\rm esc}$) parameterize unresolved ISM/feedback processes.

RT is the dominant computational expense for self-consistent, multi-source, multi-frequency RHD because it couples local radiation fields to cell/particle-scale physics everywhere.

2. Numerical Methodologies and Solvers

RHD codes employ a diverse set of approaches:

• Grid-based AMR codes: RAMSES-RT implements multi-group RT and M1 closure (Rosdahl et al., 2013); Enzo+Moray uses adaptive ray-tracing and photon-conserving schemes (Wise et al., 2010). Both codes iteratively refine grid cells in high-density regions, yielding high spatial/temporal dynamic range.
• SPH-based particle codes: TRAPHIC photon-packet-cone transport is used in GADGET and Aurora (Pawlik et al., 2016, Pawlik et al., 2015), allowing RT at the native gas resolution with cost independent of source number due to the merging of photon packets.
• Hybrid CPU-GPU codes: RAMSES-CUDATON, used in the Cosmic Dawn/CoDa project, offloads the moment-based radiative transfer and non-equilibrium hydrogen ionization to GPUs while hydrodynamics/gravity evolve on CPUs, permitting extremely large unigrid volumes (Ocvirk et al., 2015).
• Time-integration strategy: Stiff source terms arising from radiative heating/cooling and chemistry are handled with operator splitting and adaptive sub-cycling—i.e., sub-stepping chemistry/RT within each (longer) hydro time step (Whalen et al., 2010).
• Spectral treatment: Until recently, most RHD was “grey” (one/few frequency bins); modern solvers (e.g., Gasoline2+Trevr2) implement narrow-band, piece-wise power-law (PWPL) spectral reconstruction to achieve high-fidelity photoionization/heating rates, critical for CGM/ISM studies (Baumschlager et al., 2023). Piece-wise constant band approximations lead to ionization/heating errors $\gg 1$ dex in some regimes.

RT frameworks trade off between accuracy, cost, scalability, treatment of multiple sources, and ability to resolve both mean free paths and sharp ionization fronts.

3. Parameter Calibration, Convergence, and Statistical Models

Subgrid models (star formation efficiency, SN/feedback energetics, unresolved escape fractions, metal mixing) are calibrated to key observables:

• Star formation rate (SFR) functions: E.g., Aurora models tune SN feedback so that SFR functions match at $z\sim7$ above SFR$\gtrsim0.2\,M_\odot$/yr (Pawlik et al., 2016).
• Epoch of Reionization (EoR) diagnostics: $\tau_e$, CMB Thomson optical depth, reionization midpoint redshift $z_{0.5}$, and ionization duration $\Delta z$ are enforced/compared with Planck/other data (Ocvirk et al., 2015, Doussot et al., 2017, Pawlik et al., 2015).
• Convergence testing: Varying mass/spatial resolution and box size to check statistical and physical robustness. For RHD, convergence of IGM clumping, Lyman-limit system statistics, SFRD, and reionization history with increasing resolution have been systematically quantified (Pawlik et al., 2015, Pawlik et al., 2016).

Statistical models—most notably scale-dependent density-reionization bias—map expensive RadHydro outcomes to lower-resolution or larger-volume simulations:

• Battaglia et al. model: Maps the correlation between matter density $\delta_m$ and the reionization-redshift field $\delta_z$ via a scale-dependent bias $b_{zm}(k)$:

$$b_{zm}(k) = A\left[1 + (k/k_0)^{\alpha}\right]^{-1}$$

enabling fast synthesis of reionization fields that preserve the statistics of direct RHD (Battaglia et al., 2012).

4. Physical Insights and Astrophysical Results

Cosmological RHD simulations have established the following:

• Patchiness and timing of reionization: IGM is reionized “inside-out,” with ionized bubbles correlated to galaxy clustering; late/EoR timing and duration constraints are sensitive to the redshift dependence of $f_{\rm esc}$ (Doussot et al., 2017).
• Suppression of low-mass galaxy formation: Photoheating raises the IGM Jeans mass ($M_{\rm supp}\sim 2\times10^9\,M_\odot$), quenching SFRs in minihaloes and alleviating the “missing satellites” problem (Ocvirk et al., 2015).
• Feedback effects: SN feedback dominates SFR suppression across the bulk of galaxy mass function, with photoheating more important in the lowest-mass haloes. Both effects are nonlinear and coupled; omission of either yields quantitatively incorrect reionization/galaxy statistics (Pawlik et al., 2015).
• IGM/CGM structure and thermodynamics: Spectrally-accurate RT reveals extended neutral reservoirs, realistic self-shielded structures, spectral hardening, and the correct populations of UV metal ion absorbers (e.g., OVI, NV) (Cen et al., 2016, Baumschlager et al., 2023).
• Ionization-front instabilities: Nonlinear thin-shell, shadow, and angle-of-incidence instabilities at the I-fronts drive clumping, enhance UV escape fractions, and may trigger secondary star formation in protected dense clumps (Whalen et al., 2010).

Empirical sub-grid star formation models anchored to high-resolution RHD outputs enable simulating $\gtrsim (500\,\mathrm{cMpc})^3$ volumes required for statistical 21 cm predictions (Gillet et al., 2021).

5. Applications and Observational Connections

RHD outputs are used to synthesize or interpret:

• 21 cm line statistics: Full light-cone volumes at 1 Mpc resolution are produced for global signal/power spectrum forecasts (e.g., LOFAR/NENUFAR observability), incorporating spatial fluctuations in spin temperature, thermal state, and ionization (Gillet et al., 2021).
• Mock CMB secondary anisotropies: Simulated kSZ maps derive from velocity-weighted ionization structure at reionization (Battaglia et al., 2012).
• Synthetic absorption/emission lines: Metal-line and HI Lyman-alpha statistics, CGM/IGM transmission (Cen et al., 2016).
• Galaxy/halo catalogs: Star-formation histories, luminosity functions, and metallicities are traced to the RHD subgrid models and feedback-impacted formation (Ocvirk et al., 2015, Pawlik et al., 2016).
• Parameter studies: Varying $f_{\rm esc}$, source models, and feedback prescriptions assesses allowed degeneracies and matches to indirect reionization duration and symmetry constraints (Doussot et al., 2017).

Fast, bias-based mapping models provide essential support for making theoretical predictions and generating mock observations at computationally prohibitive volumes.

6. Outstanding Challenges and Future Directions

Major current challenges and open directions include:

• Multi-frequency, metal, and dust physics: Most large-scale RHD has used grey RT and primordial chemistry. Modern codes are adding detailed spectral treatment (PWPL, non-equilibrium metals, dust absorption/scattering) (Baumschlager et al., 2023).
• Resolution versus volume: Accurate treatment of both large-scale topology (mock 21 cm cubes) and small-scale clumping/radiative transfer remains computationally prohibitive. Subgrid SFR and analytic bias proxies are bridging the gap (Gillet et al., 2021, Battaglia et al., 2012).
• Reduced speed-of-light approximation (RSLA): To make explicit RT integration feasible, RSLA is widely adopted. Its limits and errors for various reionization scenarios have been quantified and are a focus of ongoing performance tests (Rosdahl et al., 2013).
• Calibration and degeneracies: The need for empirical calibration of $f_{\rm esc}$, SN feedback, and unresolved absorption fractions remains a principal uncertainty. Observational advances in metal absorption systems, 21 cm constraints, and resolved local faint galaxies are expected to progressively reduce these uncertainties (Pawlik et al., 2016, Cen et al., 2016).
• Extreme feedback and rare objects: Incorporation of AGN, X-ray binaries, soft/hard X-ray feedback, and supermassive black holes—rare but potentially dominant in some regimes—is a target for next-generation RHD (Baumschlager et al., 2023).
• Instabilities and non-linearities: I-front instabilities and turbulence in coupled RHD/chemistry/MHD frameworks are only partially resolved (Whalen et al., 2010).

A plausible implication is that the future of cosmological RHD will increasingly entail hybrid approaches—combining high-resolution subgrid modeling, large-volume bias parameterizations, and efficient RT solvers with spectrally-resolved multi-physics.

7. Summary Table: Major Codes and Methodological Approaches

Code/System	RT Method	Hydro Framework	Key Features
RAMSES-RT	M1 moments	AMR Eulerian	Multi-group, RSLA, AMR, Scalable (Rosdahl et al., 2013)
Enzo+Moray	Adaptive rays	AMR Eulerian	Photon-conserving, splitting/merging (Wise et al., 2010)
TRAPHIC (Gadget, Aurora)	Cone packets	SPH Lagrangian	Adaptive, source-agnostic RT (Pawlik et al., 2016)
RAMSES-CUDATON	M1 (GPU)	Unigrid/AMR (CPU)	Massive hybrid, pre-constrained ICs (Ocvirk et al., 2015)
Gasoline2+Trevr2	PWPL spectral	SPH Lagrangian	Narrow-band, analytic RT, CGM/ISM (Baumschlager et al., 2023)
EMMA	M1 moments	AMR	Coarse grid, empirical SFR; 21 cm focus (Gillet et al., 2021)

Key choices—hydrodynamics scheme, RT technique, spectral fidelity, feedback and subgrid calibration—dictate the scientific reach and computational tractability of each implementation, and determine the faithfulness of the simulation outputs for comparison with astrophysical data.

---

Cosmological radiation-hydrodynamic simulations are foundational tools for theoretical and observational cosmology. Their continuous evolution in methodology, physical realism, and computational design is driven both by advances in supercomputing and by the increasing stringency of multiwavelength astrophysical constraints.