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

Andrei Mesinger; Jordan Flitter; Julian B. Mu\~noz

REVIEW 2 major objections 4 minor 2 cited by

This paper claims that Lyman-alpha multiple scattering, previously requiring expensive radiative-transfer codes, can be captured by a single-parameter beta distribution, and that including it increases the 21-cm power spectrum by about 50%

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

T0 review · deepseek-v4-flash

2026-08-03 09:14 UTC pith:MZBKKWUQ

load-bearing objection Solid, transparent semi-analytic method for Lyα multiple scattering, but the headline ~50% high-redshift effect sits on an extrapolated fit in the regime the paper says dominates — don't trust the number yet. the 2 major comments →

arxiv 2601.14360 v2 pith:MZBKKWUQ submitted 2026-01-20 astro-ph.CO

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

Jordan Flitter , Julian B. Mu\~noz , Andrei Mesinger This is my paper

classification astro-ph.CO

keywords Lyman-alpha multiple scattering21-cm power spectrumcosmic dawnbeta distribution window functionWouthuysen-Field couplingsemi-numerical simulationdiffusion scaleLyman-alpha heating

verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

During the cosmic dawn, Lyman-alpha photons from the first galaxies scatter repeatedly before being absorbed, but most fast simulations pretend they travel in straight lines. The paper aims to show that the true spreading of these photons can be encoded in a single analytical window function, built from a beta distribution whose shape is set by one number: the ratio of the straight-line travel distance to a characteristic diffusion scale. Implementing this window function in a standard semi-numerical simulation, the paper finds that multiple scattering makes the 21-cm power spectrum about 50% larger at z=20, when the spin temperature is not yet locked to the gas. By z=11 and below the effect vanishes, and Lyman-alpha heating is not amplified by multiple scattering. If right, this gives a cheap way to include a physical effect that previously demanded radiative-transfer simulations.

Core claim

The central discovery is that the distance distribution of Lyman-alpha photons relative to the point where they are absorbed is universal when scaled by the diffusion scale, and is well described by a beta distribution determined by the single parameter x_em = R_SL/R_*. The resulting window function, expressed as a hypergeometric function, reduces to the straight-line window function at large scales and to a delta function at small scales, with a smooth transition in between. Because it preserves more power on small scales and decays faster on large scales than the top-hat filter, it enhances the small-scale contrast of the Lyman-alpha flux, which translates into a roughly 50% increase in th

What carries the argument

The load-bearing object is the comoving diffusion scale R_*, the distance within which a Lyman-alpha photon is likely to scatter multiple times before absorption. Its ratio to the straight-line distance R_SL defines x_em, which controls the shape of the beta distribution via empirical fits for the mean and variance. The Fourier transform of that beta distribution gives the multiple-scattering (MS) window function, a generalized hypergeometric function, which replaces the straight-line top-hat filter when convolving the emissivity field to compute the Lyman-alpha flux.

Load-bearing premise

That the beta-distribution parameters extrapolated to x_em<0.2, the regime that contributes most to the Lyman-alpha flux at z=20, are accurate enough for the quantitative 50% power-spectrum result.

What would settle it

Run a full radiative-transfer simulation at z=20 in a comparable 300 Mpc box that includes non-linear peculiar velocities and finite temperature, compute the 21-cm power spectrum with and without the straight-line approximation, and compare the ratio to the predicted ~50% increase. Alternatively, measure the radial distribution of absorbed photons at x_em<0.2 directly with high-statistics Monte Carlo and check whether a beta distribution with the paper's extrapolated parameters fits the data.

Watch this falsifier — get emailed when new claim-graph text bears on it.

If this is right

Semi-numerical 21-cm simulations can now include Lyman-alpha multiple scattering at virtually no extra computational cost, enabling large parameter-space exploration for cosmic dawn.
Parameter constraints from high-redshift (z~20) 21-cm observations would shift by the ~50% power-spectrum change if multiple scattering is ignored, biasing inferred star formation and X-ray properties.
At redshifts below about 15, the straight-line approximation remains valid for the 21-cm power spectrum, so existing low-redshift measurements do not need revision.
The combined effect of multiple scattering and Lyman-alpha heating is negligible, so the two processes can be treated independently in simulations.
The window function's redshift dependence means that using the straight-line filter across all redshifts misestimates the early-time Lyman-alpha flux.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

The single-parameter universality suggests the same window function may hold under modified cosmologies or different astrophysical models, because the diffusion scale absorbs most of the redshift and ionization dependence; the paper only offers a cautious prediction here.
The paper's disagreement with earlier work on Lyman-alpha heating implies the equilibrium heating calculation is sensitive to implementation details; a controlled code comparison would clarify whether the discrepancy comes from the window function or the heating rates.
If the x_em<0.2 extrapolation is wrong, the magnitude of the 50% effect at z=20 would shift, but the qualitative conclusion that multiple scattering matters only at high redshifts would likely survive.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Desk Editor's Note

Solid, transparent semi-analytic method for Lyα multiple scattering, but the headline ~50% high-redshift effect sits on an extrapolated fit in the regime the paper says dominates — don't trust the number yet.

read the letter

Jordan, quick note on Flitter et al. 2601.14360. The useful news first: the paper gives the 21-cm community a genuine tool rather than another black box. The diffusion-scale argument from Loeb-Rybicki is clean, the single-parameter x_em compression of the MC radial distributions is plausible, and the beta-distribution window function with finite-shell correction is worked out carefully. The code is public, the implementation details in Appendix D (series + asymptotic matching) are exactly what someone needs to reproduce it. The limit x_em→∞ returning the straight-line window function is a nice sanity check. On those grounds the paper deserves a serious referee.

The soft spot is exactly where the reader and stress-test point: the ~50% power-spectrum enhancement at z=20 is computed in a regime where the fit is an extrapolation. Equations (29)-(30) are explicitly power-law extrapolations below x_em=0.2, and Sec. 5 says small filter radii dominate the J_alpha integral at high redshift. With R_*≈31 Mpc at z=20, the smallest 21cmFAST shell at R_i≈0.93 Mpc corresponds to x_em≈0.03. So the dominant contribution to the flux at the redshift where the effect matters comes from x_em values that were never validated against SPαRTA output. That doesn't automatically sink the qualitative conclusion — 'only at high z, no heating boost' is probably right — but the 50% number is a prediction of an extrapolated function, not a measurement. I'd want to see additional MC samples at x_em<0.2 or a cross-check against a full radiative-transfer run before quoting that number.

Other soft spots are minor and mostly acknowledged: N_photons=10^3 is low (they claim convergence, no error bars shown on the 50%), density inhomogeneities are set to zero, and the peculiar-velocity treatment is linear, which is why they diverge from Semelin et al. and Mittal et al. Those are limitations, not misrepresentations. The citation pattern is fair; the RFB21 six-parameter fit is properly distinguished.

Verdict: the formalism is a contribution and the paper is honest about its assumptions. The unvalidated x_em<0.2 extrapolation is a load-bearing flaw for the headline number but not for the method. I'd accept it for peer review without a second thought, with a request for the small-x_em validation. Would I cite it? Likely yes, for the window function and the implementation, and I'd flag the 50% as provisional. Bring it to reading group? Yes — the extrapolation-vs-dominance tension is a good discussion.

Referee Report

2 major / 4 minor

Summary. The paper develops a semi-analytical treatment of Lyα multiple scattering (MS) for use in semi-numerical 21-cm simulations. The authors introduce SPαRTA, a fast Monte Carlo code that tracks photon trajectories in an expanding IGM, and show that the radial distribution of absorbed Lyα photons is well described by a beta distribution controlled by a single parameter x_em, the ratio of the straight-line distance to the comoving diffusion scale R_*. They derive an analytic window function (Eq. 32) for the MS filtering of the emissivity field and implement it in 21cmFAST. Their central quantitative result is that including MS changes the 21-cm power spectrum by about 50% at z=20, when Lyα coupling is inefficient, while being negligible at lower redshifts and for Lyα heating. The paper also reports a tension with earlier work by RFB21 regarding Lyα heating.

Significance. If the result holds, the paper provides a computationally efficient way to include a physical effect that has previously required expensive radiative-transfer simulations. The analytic window-function derivation, the public release of SPαRTA, and the explicit implementation in 21cmFAST are concrete strengths. The calibrated beta-distribution fit and the asymptotic matching to the straight-line limit are well executed. The claimed ~50% effect at z=20, if robust, is observationally relevant for upcoming 21-cm experiments. However, as detailed below, the quantitative headline result relies on an extrapolated regime of the fit, and the lack of Monte Carlo error estimates weakens the current support for that number.

major comments (2)

[Sec. 4.3 (Eqs. 29–30) and Sec. 5 (Fig. 8)] The central ~50% power-spectrum claim at z=20 is driven by small filter radii, as Sec. 5 states that contributions from small radii dominate the J_α integral. The smallest shells shown in Fig. 8 (R_i≈0.93 Mpc) correspond to x_em≈0.03, well below the x_em≥0.2 range where the fits of Eqs. (29)–(30) are validated. The text explicitly calls the x_em≤0.2 branch a power-law extrapolation. Since this unvalidated regime dominates the z=20 signal, the quantitative result is not yet supported. Please either run SPαRTA at x_em<0.2 and validate the extrapolation, or demonstrate that the ~50% effect is robust to plausible alternative extrapolations of μ(x_em) and η(x_em).
[Sec. 3 (N_photons=1000) and Figs. 3–6, 9–11] No Monte Carlo statistical uncertainties are shown for the fitted μ and η, nor in the propagated power-spectrum differences. The paper states that results are insensitive to N_photons and L_cell and claims convergence, but no quantitative evidence is provided. Given that the headline result is a ~50% enhancement and that Fig. 9 presents it at all k, the absence of any error bars or convergence plot makes the claimed amplitude difficult to assess. Please add error bars to the beta-distribution fits and to the resulting Δ²_21 differences, and show the convergence test for N_photons and L_cell.

minor comments (4)

[Eqs. (29)–(30)] The variable ζ_em is used without definition in the equation block; please define it (ζ_em≡log₁₀ x_em) in the text or immediately below the equations.
[Eq. (4)] The approximation f_v≈δ_D is stated and deferred to Sec. 4 for justification. It would improve readability to add a one-sentence justification here, e.g., that v_rel/c≪1 for the relevant scales.
[Fig. 8] The phrase 'front side of the SFRD coeval box' is unclear; specify whether this is a slice, a projection, or a face of the 3D box.
[Appendix D.1] Eq. (D.3) uses the notation '1F1' and 'ix' in a way that could be misread as the imaginary unit; since x is real, this is fine, but consider writing ±ix explicitly or using z = ix to avoid confusion.

Circularity Check

0 steps flagged

No circularity: the beta-distribution fit is calibrated to independent SPαRTA Monte Carlo output, and the 21-cm power-spectrum result is a downstream application, not a re-fit or self-referential prediction.

full rationale

The derivation chain is not circular. The diffusion scale R_*(z_abs) is obtained independently from Loeb & Rybicki (1999) via Eq. (24), and the beta-distribution parameters µ(x_em) and η(x_em) are fit to SPαRTA Monte Carlo radial-distribution data (Sec. 4.2–4.3), not to the 21-cm power spectrum. The MS window function, Eq. (32), is a direct Fourier/hypergeometric transform of that fitted radial distribution, and its insertion into 21cmFAST is a forward propagation: the ~50% power-spectrum enhancement at z=20 emerges from filtering the SFRD field, with no equation reducing that enhancement to the fitted constants by construction. The acknowledged power-law extrapolation for x_em≤0.2 (Sec. 4.3) combined with the statement that small filter radii dominate J_α at z=20 is a genuine validity/correctness risk, but it is not circularity: the prediction is not used to adjust the fit, and the fit is fixed by data at x_em≥0.2 plus an explicit extrapolation assumption. The paper's self-citations (e.g., EOS2026 parameters from Breitman et al. 2026, Lyα-heating implementation from Sarkar et al. 2022) are not load-bearing for the MS window function claim: the EOS2026 parameters are anchored to external JWST/reionization data, and the heating implementation supports the subsidiary conclusion that Lyα heating remains negligible. No uniqueness theorem, no ansatz-by-citation, and no renamed known result is invoked to force the central claim. The honest non-finding is therefore score 0, with the extrapolation concern noted as a limitation rather than circularity.

Axiom & Free-Parameter Ledger

5 free parameters · 6 axioms · 0 invented entities

The central claim is not derived from first principles alone. It rests on a Monte Carlo code with several chosen numerical parameters, an assumed beta functional form, an extrapolated small-x_em fit, and neglect of density and non-linear-velocity effects. These are the price of admission for a fast semi-analytic implementation, but they limit the certainty of the quantitative 50% result.

free parameters (5)

Beta-fit parameters μ(x_em), η(x_em) = Piecewise coefficients in Eqs. (29)-(30)
Fit to SPαRTA output; fully determines the multiple-scattering window function via Eq. (32).
First-step redshift interval Δz/(1+z_abs) = 2×10^-4
Chosen to start the Monte Carlo outside the Voigt core; sets the initial photon frequency.
Virtual cell size L_cell = 200 kpc
Numerical resolution choice for step size; convergence asserted but not shown.
Photon number N_photons = 1000
Low raw photon count; intermediate points are reused as effective sources. Convergence asserted, not demonstrated in a figure.
IGM temperature T_k and ionization fraction x_HI in SPαRTA = T_k = 10^4 K, x_HI = 1 (fiducial)
Physical inputs to the Monte Carlo; the paper argues the radial distributions collapse to x_em alone, but spatial fluctuations of T_k and x_HI are not modeled.

axioms (6)

domain assumption The spatial photon distribution f is homogeneous, isotropic, and independent of the absorption point (Eq. 5).
Enables a single Fourier-space window function; the paper notes it breaks down during patchy reionization and argues it has little effect on brightness temperature.
ad hoc to paper The radial distributions are exactly beta distributions on y∈[0,1] (Eq. 26).
Chosen for finite support and single-peak shape; no physical derivation of the beta form is given.
domain assumption IGM density inhomogeneities are neglected (δ=0) in SPαRTA (Appendix C.1).
A joint draw of δ and v_parallel from one past sample yields a non-positive-definite covariance matrix; justified by citing Mittal et al. 2024.
domain assumption Linear perturbation theory with conditional Gaussian velocities from a single past sample is sufficient for peculiar velocities (Appendices B-C).
Ignores non-linear velocity enhancements that Semelin et al. 2023 and Mittal et al. 2024 argue are important; the paper explicitly defers this to future work.
domain assumption Matter domination, Voigt-wing cross-section, and frozen-redshift approximations in the diffusion-scale derivation (Eqs. 22-24).
Used to derive the analytic R_* that defines x_em; valid for the z~5-35 regime but approximate.
domain assumption Doppler broadening of the apparent SED is approximated as a delta function in relative velocity (Eq. 4).
Assumes v_rel/c << 1 and a flat SED; the paper says this is justified for flat SEDs.

pith-pipeline@v1.3.0-alltime-deepseek · 34098 in / 13974 out tokens · 152315 ms · 2026-08-03T09:14:50.814857+00:00 · methodology

0 comments

read the original abstract

A crucial physical quantity in determining the 21-cm signal during cosmic dawn is the inhomogeneous background of Ly$\alpha$ 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$\alpha$ and Ly$\beta$ are assumed to travel in straight lines until they redshift into the Ly$\alpha$ resonance. However, due to the damping wing in the Ly$\alpha$ 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$\alpha$ 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$\alpha$ multiple scatterings in the public, semi-numerical simulation 21cmFAST. We employ Monte Carlo simulations to study the trajectories of Ly$\alpha$ photons on different scales. We find that the distance distributions of Ly$\alpha$ 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$\alpha$ multiple scattering does not enhance Ly$\alpha$ heating, and that the combined effect is negligible, especially under realistic X-ray heating scenarios.

Figures

Figures reproduced from arXiv: 2601.14360 by Andrei Mesinger, Jordan Flitter, Julian B. Mu\~noz.

**Figure 1.** Figure 1: Top row: trajectories of five photons, as was simulated by SPαRTA. All photons are absorbed at the center once their apparent frequency reaches Lyα at zabs = 10. It is assumed that the photons travel in a homogeneous IGM with Tk = 104 K and xHI = 1 (though the gas is allowed to move with some peculiar velocity with respect to the comoving frame). For reference, we show by the black line the trajectory of a… view at source ↗

**Figure 2.** Figure 2: Left panel: Radial distance of the photons from the absorption point (at zabs = 10) as a function of redshift. Right panel: Apparent frequency shift of the photons from Lyα as a function of redshift. The colors of the curves here match the colors of the trajectories in [PITH_FULL_IMAGE:figures/full_fig_p007_2.png] view at source ↗

**Figure 3.** Figure 3: Radial distributions fMS (y;zabs) of the normalized distance variable, y ≡ r/RSL, for different values of xem ≡ RSL (zabs,zem) /R∗ (zabs) (for better interpretation of that quantity, we report also the corresponding ∆z = zem − zabs). All curves shown here were obtained by fitting the numerical distributions from SPαRTA to a beta function. Solid (dashed) curves correspond to having turned on (off) all the… view at source ↗

**Figure 4.** Figure 4: Left panel: distribution of the relative peculiar velocity between the absorber at zabs = 10 and the emitter (whose redshift is implied by the xem value), as extracted from the SPαRTA simulation. All curves are cut at the highest and lowest values obtained in the simulation. We also show, by the dashed black curve, the Gaussian distribution that corresponds to |v|/c at z = 0 (computed in Newtonian gauge), … view at source ↗

**Figure 5.** Figure 5: Similarly as [PITH_FULL_IMAGE:figures/full_fig_p010_5.png] view at source ↗

**Figure 7.** Figure 7: The MS window function, as given by Eq. (32) and the fit of Eqs. (28)-(30), for different xem values. For comparison, we show the SL window function, Eq. (13), by the black dashed curve. For better visualization, the inset shows how the window functions behave at high kR values. A few remarks on the fit of Eqs. (29)-(30). Firstly, the fit for µ (xem) is much more accurate compared to η (xem), and matches t… view at source ↗

**Figure 8.** Figure 8: Comparison between SL filtering and MS filtering. Top row shows (from left to right) the front side of the unfiltered SFRD coeval box at z = 20 and the power spectrum of the box. Middle row shows (from left to right) the filtered SFRD box with the SL (MS) window function, Eq. (11) and Eq. (13) (Eq. 32), and a comparison of the power spectra of the two filtered boxes, for Ri = 0.93 Mpc, Ro = 1.09 Mpc. Botto… view at source ↗

**Figure 9.** Figure 9: Comparison between SL filtering and MS filtering on the relevant physical fields at z = 20. Top row shows (from left to right) the front side of the coeval box of Jα (z) /J0 (z) at z = 20 when the SL filter is applied on the SFRD field (see [PITH_FULL_IMAGE:figures/full_fig_p014_9.png] view at source ↗

**Figure 10.** Figure 10: Similarly as [PITH_FULL_IMAGE:figures/full_fig_p015_10.png] view at source ↗

**Figure 11.** Figure 11: Left panel: evolution of the mean Tk as a function of redshift, for log10 (LX/SFR) = 35 and two different star formation efficiencies (in our fiducial model, f (II) ∗,fid = 1%, f (III) ∗,fid = 0.03%). The black curve corresponds to a scenario with no Lyα heating, while the blue (orange) curves account for Lyα heating, with SL (MS) filtering. Shaded regions around the mean value have width of the RMS of th… view at source ↗

discussion (0)

Forward citations

Cited by 2 Pith papers

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score.

When galaxies burst II. Implications of enhanced burstiness for the 21-cm Cosmic Dawn signal
astro-ph.CO 2026-07 conditional novelty 6.0

Bursty star formation adds a shot-noise-like term to the Cosmic Dawn 21-cm power spectrum, boosting it by factors of a few at the Wouthuysen-Field trough and strongly on small scales at high redshift.
Implications of \textit{SARAS3} data for Coulomb-like interacting dark matter
astro-ph.CO 2026-05 unverdicted novelty 5.0

SARAS3 data sets T21 ≳ -277.6 mK at z=23.6 (3σ) in IDM models and finds no preference for IDM over CDM (Bayes factor ≈1.7).

Reference graph

Works this paper leans on

3 extracted references · 1 linked inside Pith · cited by 2 Pith papers

[1]

Abdurashidova, Z. et al. 2022, Astrophys. J., 924, 51 Adi, T., Flitter, J., & Kovetz, E. D. 2025, Phys. Rev. D, 111, 043515 Aghanim, N. et al. 2020, Astron. Astrophys., 641, A6, [Erratum: As- tron.Astrophys. 652, C4 (2021)] Bale, S. D., Bassett, N., Burns, J. O., et al. 2023, arXiv e-prints, arXiv:2301.10345 Barkana, R. 2018, Nature, 555, 71 Barkana, R. &...

Pith/arXiv arXiv 2022
[2]

WA,B (kr),(C.3) where ˜Tδ (k,z )≡T δ (k,z ) while ˜Tvi (k,z )≡T θ (k,z )/( √ 3k), and the window functionsW A,B (kr) are listed in Table C.1. From that table we can see that two perpendicular components of the pe- culiar velocity,v⊥ 1 andv⊥ 2 , are completely uncorrelated with any other field, and thus they depend only on their own past samples. However, ...

2011
[3]

In other words, at every step in the simulation, the size of the joint covariance matrix needs to be increased until it is positive-definite

(or vice versa). In other words, at every step in the simulation, the size of the joint covariance matrix needs to be increased until it is positive-definite. While this is technically possible, we defer the implementation of this solution to future work. In the meantime, in order to overcome this mathematical ob- stacle, in this work we have completely i...

2024

[1] [1]

Abdurashidova, Z. et al. 2022, Astrophys. J., 924, 51 Adi, T., Flitter, J., & Kovetz, E. D. 2025, Phys. Rev. D, 111, 043515 Aghanim, N. et al. 2020, Astron. Astrophys., 641, A6, [Erratum: As- tron.Astrophys. 652, C4 (2021)] Bale, S. D., Bassett, N., Burns, J. O., et al. 2023, arXiv e-prints, arXiv:2301.10345 Barkana, R. 2018, Nature, 555, 71 Barkana, R. &...

Pith/arXiv arXiv 2022

[2] [2]

WA,B (kr),(C.3) where ˜Tδ (k,z )≡T δ (k,z ) while ˜Tvi (k,z )≡T θ (k,z )/( √ 3k), and the window functionsW A,B (kr) are listed in Table C.1. From that table we can see that two perpendicular components of the pe- culiar velocity,v⊥ 1 andv⊥ 2 , are completely uncorrelated with any other field, and thus they depend only on their own past samples. However, ...

2011

[3] [3]

In other words, at every step in the simulation, the size of the joint covariance matrix needs to be increased until it is positive-definite

(or vice versa). In other words, at every step in the simulation, the size of the joint covariance matrix needs to be increased until it is positive-definite. While this is technically possible, we defer the implementation of this solution to future work. In the meantime, in order to overcome this mathematical ob- stacle, in this work we have completely i...

2024