Semi-analytical approach to Ly$α$ multiple-scattering in 21-cm signal simulations

A crucial physical quantity in determining the 21-cm signal during cosmic dawn is the inhomogeneous background of Ly$α$ photons originating from the first galaxies. As these photons travel through the intergalactic medium, their scattering cross-section is often approximated as a delta function at resonance due to computational cost. That is, photons with emitted wavelengths between Ly$α$ and Ly$β$ are assumed to travel in straight lines until they redshift into the Ly$α$ resonance. However, due to the damping wing in the Ly$α$ cross-section, this approximation fails as the frequency of the photon approaches the resonant frequency, resulting in multiple scatterings events that could be separated by non-negligible distances. Some previous works studied this effect of Ly$α$ multiple scattering by running computationally heavy radiative-transfer simulations. However, robustly interpreting the cosmic 21cm signal requires exploring a large parameter space of astrophysical uncertainties, motivating more computationally-efficient approaches. Here we incorporate Ly$α$ multiple scatterings in the public, semi-numerical simulation 21cmFAST. We employ Monte Carlo simulations to study the trajectories of Ly$α$ photons on different scales. We find that the distance distributions of Ly$α$ photons with respect to the absorption point can be modeled as analytical functions that are governed by a single parameter. Upon implementing the distance distributions in 21cmFAST, we find that the multiple scattering effect is important (about 50% difference in the 21-cm power spectrum) only at high redshifts before the spin temperature is fully coupled to the kinetic temperature. Furthermore, we find that Ly$α$ multiple scattering does not enhance Ly$α$ heating, and that the combined effect is negligible, especially under realistic X-ray heating scenarios.

💡 Research Summary

This paper presents a computationally efficient, semi‑analytical method to incorporate Ly α multiple‑scattering (MS) effects into the widely used semi‑numerical 21 cm simulation code 21cmFAST. The motivation stems from the fact that the Ly α background, which drives the Wouthuysen‑Field coupling (through the coupling coefficient xα), is usually modeled with a delta‑function approximation of the Ly α cross‑section. In that “straight‑line” (SL) approximation, photons emitted between Ly β and Ly α travel unimpeded until cosmological redshift brings them exactly to resonance, at which point a single scattering occurs. However, the damping wing of the Ly α cross‑section gives a non‑zero probability of scattering before exact resonance, allowing photons to undergo several scatterings separated by finite comoving distances. This effect can modify the spatial distribution of Ly α flux Jα and consequently the 21 cm brightness temperature fluctuations, especially at high redshifts (z ≈ 20–30) where the spin temperature Ts is not yet fully coupled to the kinetic temperature Tk.

The authors develop a public Monte‑Carlo ray‑tracing code called SPαRTA (Speedy Ly α Ray Tracing Algorithm). SPαRTA does not rely on a spatial grid; instead it uses analytical results from linear perturbation theory to emulate peculiar‑velocity effects and includes finite IGM temperature. The code runs in seconds to minutes, producing photon trajectories for a wide range of IGM conditions. By simulating many photon histories, the authors find that the probability distribution of the comoving distance r between photon emission and its final absorption point follows a beta distribution:

 P(y) ∝ y^{α−1}(1−y)^{β−1}, y ≡ r/RSL,

where RSL(z, z′) is the straight‑line comoving distance corresponding to the redshift interval between emission (z′) and absorption (z). Remarkably, only a single shape parameter β varies significantly with IGM temperature and the initial photon frequency; the other shape parameter is essentially fixed. In the limit β → ∞ the beta distribution collapses to a delta function, reproducing the SL approximation.

Armed with this analytic form, the authors construct a “window function” W(k; z, z′), the Fourier transform of the radial distribution f(r; z, z′). For an isotropic, purely radial distribution the window function reduces to an integral over y of the beta distribution multiplied by sin(k RSL y)/(k RSL y). This window function can be multiplied in Fourier space with the unfiltered emissivity field ε∗(x, z′) (which encodes the star‑formation rate density and spectral energy distribution of sources) to obtain the apparent emissivity at the absorption redshift. The apparent Ly α flux Jα is then computed by integrating over source redshifts, exactly as in the standard 21cmFAST formalism, but now with the MS‑modified kernel.

The authors implement this kernel in 21cmFAST and compare results with the traditional SL treatment. They find that:

1. At high redshifts (z ≈ 25) where xα ≈ 1, the inclusion of MS raises the mean Ly α flux by ~20 % and boosts the 21 cm power spectrum P21(k) by up to ~50 % at k ≈ 0.1 Mpc⁻¹. This reflects the fact that spatial fluctuations in Jα dominate the signal when the coupling is weak.

2. Once the spin temperature fully couples (xα ≫ 1, typically z ≲ 15), the Ly α flux saturates and the MS‑induced modifications become negligible; the power spectrum differences drop below a few percent.

3. The additional heating from Ly α photons (Ly α heating) is essentially unchanged by MS because the total number of scatterings and the average energy transferred per scattering remain the same. Even when the X‑ray heating efficiency is set to very low values (well below current HERA or JWST constraints), the combined effect of MS and Ly α heating alters the 21 cm signal by less than ~5 %.

These findings have two practical implications. First, the semi‑analytical beta‑distribution kernel provides a physically motivated yet computationally cheap way to include Ly α MS in large‑scale parameter‑space explorations (e.g., MCMC, neural‑network emulators) without resorting to full radiative‑transfer simulations, which are orders of magnitude more expensive. Second, for upcoming 21 cm experiments (HERA, SKA), the MS effect only needs to be accounted for when interpreting the high‑redshift (cosmic‑dawn) power spectrum; for the later reionization era the standard SL approximation remains sufficient.

In summary, the paper delivers a robust, open‑source tool (SPαRTA) and an analytic prescription (beta‑distribution window function) that together enable efficient inclusion of Ly α multiple‑scattering in semi‑numerical 21 cm simulations. The approach bridges the gap between fully detailed radiative‑transfer calculations and the speed required for modern Bayesian inference, while clarifying that, under realistic astrophysical scenarios, the impact of MS on the observable 21 cm signal is modest and confined to the earliest epochs of cosmic dawn.