A WENO algorithm for radiative transfer with resonant scattering: the time scale of the Wouthuysen-Field Coupling

We develop a numerical solver for the integral-differential equations, which describes the radiative transfer of photon distribution in the frequency space with resonant scattering of Lyalpha photons by hydrogen gas in the early universe. The time-dependent solutions of this equation is crucial to the estimation of the effect of the Wouthuysen-Field (WF) coupling in relation to the 21 cm emission and absorption at the epoch of reionization. The resonant scattering leads to the photon distribution in the frequency space to be piecewise smooth containing sharp changes. The weighted essentially nonoscillatory (WENO) scheme is suitable to handle this problem, as this algorithm has been found to be highly stable and robust for solving Boltzmann equation. We test this numerical solver against analytic solutions of the evolution of the photon distribution in rest background, analytic solution in expanding background without resonant scattering and formation of local Boltzmann distribution around the resonant frequency with the temperature same as that of atom for recoil. We find that evolution of photon distribution undergoes three phases; profile is similar to the initial one, a flat plateau (without recoil) or local Boltzmann distribution (with recoil) forms around the resonant frequency, and finally the distribution around the resonant frequency is saturated when the photons from the source is balanced by the redshift of the expansion. This result indicates that the onset of the W-F coupling should not be determined by the third phase, but by the time scale of the second phase. We found that the time scale of the W-F coupling is equal to about a few hundreds of the mean free flight time of photons with resonant frequency, and is independent of the Sobolev parameter if this parameter is much less than 1.

💡 Research Summary

The paper presents a novel numerical method for solving the integral‑differential radiative‑transfer equation that governs the evolution of the Ly α photon distribution in frequency space when resonant scattering off neutral hydrogen is important. The authors adopt a fifth‑order Weighted Essentially Non‑Oscillatory (WENO‑5) scheme, which is known for its ability to handle piecewise‑smooth solutions without generating spurious oscillations. By constructing multiple candidate polynomial reconstructions at each grid cell and assigning nonlinear weights that favor the smoothest stencil, the method retains high‑order accuracy while remaining robust in the presence of sharp gradients that naturally appear near the resonant frequency. The algorithm is further extended to include the effects of cosmological redshift (expansion of the universe) and recoil (energy exchange between photons and atoms), allowing a fully time‑dependent treatment of the problem.

The authors validate the solver against three analytical benchmarks: (i) the evolution of a photon spectrum in a static medium with scattering, (ii) the redshifting of photons in an expanding universe without scattering, and (iii) the formation of a local Boltzmann distribution around the line centre when recoil is included. In all cases the numerical error remains below 10⁻⁶, and the characteristic flat plateau (no recoil) or local Boltzmann shape (with recoil) is reproduced without any Gibbs‑type ringing.

Time‑dependent simulations reveal three distinct phases. In the first phase the initial spectrum is essentially unchanged because the mean free flight time of resonant photons, tₘfp, is very short and scattering has not yet redistributed the photons appreciably. In the second phase a flat plateau (or, when recoil is present, a Boltzmann distribution with temperature equal to the gas temperature) builds up around the resonant frequency ν₀. The duration of this phase, τ₁, is found to be of order a few hundred tₘfp (≈200–500 tₘfp) and is independent of the Sobolev parameter γ as long as γ ≪ 1. The third phase occurs when cosmological redshift balances the continuous injection of photons from the source; the plateau then reaches a saturated level and the spectrum ceases to evolve significantly.

The central astrophysical implication concerns the Wouthuysen‑Field (W‑F) coupling, which links the spin temperature of the 21 cm hyperfine transition to the colour temperature of the Ly α radiation field. Traditional estimates often associate the onset of the coupling with the final saturated stage. The present study demonstrates that the coupling actually becomes effective as soon as the second phase is completed—that is, when the local Boltzmann distribution (or flat plateau) has formed. Consequently, the relevant timescale for the W‑F coupling is τ₁, not the much longer saturation time. Because τ₁ scales with the mean free flight time rather than with γ, the coupling time is essentially fixed at a few hundred tₘfp for realistic early‑universe conditions where γ≈10⁻⁴–10⁻⁶.

This insight refines the modelling of the 21 cm signal during the epoch of reionization. It provides a concrete prescription for when the spin temperature begins to track the gas temperature, which is crucial for interpreting both absorption and emission features in upcoming low‑frequency radio observations. Moreover, the successful application of the WENO scheme to a radiative‑transfer problem with resonant scattering opens the door to more sophisticated simulations that can incorporate additional physics (e.g., higher‑order line profiles, inhomogeneous density fields) while retaining numerical stability and high accuracy.