How X-rays heat the IGM in different 21-cm simulation codes: a comparison between Licorice and Beorn

Abstract: Any interpretation of the 21-cm signal of neutral hydrogen using Bayesian inference methods can only be as accurate as the underlying simulation code used to model the state of the intergalactic medium (IGM). 3D radiative transfer (RT) simulation codes may capture complex physics, but are computationally expensive and, therefore, faster, more approximate codes have been developed. To improve our understanding of the convergence of simulation codes in the 21-cm science community, we present a comparison of the X-ray heating of the IGM modelled in Licorice, a 3D RT simulation code, and Beorn, a 1D RT code. We use Beorn to process sources extracted from Licorice simulations, using the same physics of the sources, in order to obtain two versions of the temperature of the IGM heated by X-rays. We observe a good agreement between the luminosity fields, mean temperatures, and global 21-cm signal of the two setups, but discrepancies in the distribution of temperature and 21-cm signal, which result in a $\sim 30\%$ difference in the 21-cm power spectrum. We attempt to isolate the approximations that lead to these differences and find that common approximations used in 1D RT codes produce effects of that magnitude. Using an emulator of the Licorice power spectra in an MCMC pipeline, we translate these differences between power spectra into differences between posterior distributions over the astrophysical parameters. We observe a typical bias between 1D posteriors of $\gtrsim 1 σ$ (with a noise level corresponding to 100h of SKA observations).

• The paper compares identical source populations in Licorice and Beorn, finding mean-temperature differences below 5% in idealised tests and about 10% with native physics.
• The simulations produce narrower temperature distributions in Beorn, distinct heating topology, and 21-cm power-spectrum differences averaging 20% realistically and 11% under matched assumptions.
• These discrepancies can bias inferred astrophysical parameters by typically more than 1σ, with mean-density optical-depth calculations alone generating substantial errors relevant to SKA analyses.

Motivation and scope

Interpretation of 21-cm measurements through Bayesian inference is only as reliable as the simulation codes used to model the intergalactic medium (IGM). The community employs a spectrum of codes — semi-numerical tools such as 21cmFAST, 1D radiative transfer (RT) codes that paint pre-computed profiles around sources (BEARS, Grizzly, Beorn), and full 3D RT codes (C2-ray, EMMA, Licorice) — and prior comparison studies have found power-spectrum differences of order $\sim$20% between 3D RT and faster schemes. These earlier comparisons, however, focused on the ionisation field under the assumption of a saturated spin temperature, leaving X-ray heating unexamined. This paper fills that gap by comparing the X-ray heating of the IGM in Licorice, a coupled SPH–Monte Carlo ray-tracing hydro-radiative code, and Beorn, a 1D RT code that paints pre-computed temperature profiles around sources (2608.14423).

The practical stakes are considerable. The Loreli II database of nearly 10,000 Licorice simulations at $256^3$ resolution consumed several million CPU-hours, yet its sampling density is already insufficient for noise levels below roughly 100 hours of SKA observations. If cheaper codes reproduce 3D RT results accurately, expensive databases are hard to justify; if they do not, systematic modelling errors will propagate directly into posteriors over astrophysical parameters.

Matching source populations across codes

To make the comparison meaningful, the authors feed identical source populations to both codes. Since Licorice does not produce halo catalogues and its source model carries no explicit halo-mass dependence, each non-zero pixel of the UV luminosity cube is treated as a single source. Luminosity histories are constructed by tracking every source backward in time from the final redshift, handling cell crossings and mergers heuristically, and sorting sources into $N = 1000$ bins by their final luminosity. Mean histories per bin are converted into star-formation histories via the standard X-ray-to-SFR proportionality, so that Beorn's sources match Licorice's X-ray luminosities.

Two technical points deserve emphasis. First, an initial implementation in which the number of sources per bin varied with redshift caused Beorn to violate energy conservation at the profile-painting stage, producing spurious 10–30% differences in mean temperature; including zero-luminosity "unborn" sources in the binning restores energy conservation. Second, the binning procedure narrows the luminosity distribution relative to Licorice, suppressing both very bright and very dim sources — an unavoidable approximation shared broadly by 1D RT codes.

Both setups use soft X-rays only (power-law index 1.6 between 100 eV and 2 keV), with hard X-rays disabled in Licorice, and no UV emission, isolating the heating physics. Total X-ray luminosities agree at the $\sim$1% level, and sources occupy identical positions.

Structural differences between the codes

Beyond the binning approximation, three physical differences distinguish the realistic versions of the codes:

• Homogeneous density: Beorn computes optical depths assuming the mean density of the Universe, whereas Lagrangian Licorice uses local densities; heat transport by moving SPH particles is also absent in Beorn.
• Heating fraction: Licorice evaluates $f_{heat}$ using the local ionised fraction following the Shull & van Steenberg fit; Beorn uses the box-averaged electron fraction, since profiles are pre-computed independently of their environment.
• Flux overlap: 1D codes handle the overlap of heated regions around neighbouring sources imperfectly, which can significantly bias $f_{heat}$ where heating regions intersect.

These approximations are not idiosyncratic to Beorn: semi-numerical codes such as 21cmFAST likewise assume homogeneous density in optical-depth calculations, so the conclusions bear on families of approximations rather than one code alone.

Results: temperatures, fields, and power spectra

The authors run two experiments per model: a "realistic" setup with each code's native physics, and an "idealised" setup forcing constant $f_{heat} \approx 0.11$, mean-density optical depths, and volume-averaged temperatures in both codes. Five parameter sets varying $f_X$, $\tau_{SF}$, and $M_{min}$ are considered, focusing on relatively high $256^3$0 to make heating differences visible.

Aggregate quantities agree well. In the idealised setup, mean temperatures differ by less than $256^3$15%, validating the pipeline and energy conservation. In the realistic setup the agreement degrades modestly to $256^3$210%, and the global 21-cm signal shows larger discrepancies than the mean temperatures — evidence that spatial correlations between the temperature field, the Wouthuysen–Field coupling, and density differ between the codes even when bulk energetics match.

Field-level differences persist. Temperature probability distributions in Beorn are narrower, with fewer very hot and very cold pixels. Topological analysis via the first two Betti numbers shows consistent excess of loops ($256^3$3) in Licorice, indicating genuinely different heating morphology. Notably, these differences survive even in the idealised setup, implying that the identified approximations do not fully account for the discrepancy; the residual is attributed to the source-binning procedure or unidentified effects.

Power spectra disagree at the level that matters for inference. The mean absolute difference in the 21-cm power spectrum is 20% in the realistic setup and 11% in the idealised one, largest at $256^3$4–12 near the absorption peak of the global signal, where temperature fluctuations couple non-trivially to Lyman-$256^3$5 coupling and density. These differences exceed cosmic variance on large scales and thermal noise on moderate scales for 100 hours of SKA observation, reaching an average of $256^3$6 for the fiducial model. Discrepancies decrease toward lower redshift as heating saturates.

To isolate causes, the authors re-run Licorice with individual Beorn-like approximations: homogeneous density throughout RT, mean-density optical depth only, globally averaged $256^3$7, and all combined. Each approximation alone produces differences of a few tens of percent. Strikingly, merely computing the optical depth with the mean density — while solving the temperature equation with the true local density — already significantly biases the power spectrum. This identifies a specific, widely-shared approximation as individually sufficient to cause large code-to-code scatter.

Consequences for parameter inference

Using Loremu II, an emulator of Licorice power spectra, within an MCMC pipeline with a Gaussian likelihood including thermal noise (100 h SKA), cosmic variance, and emulator error, the authors infer posteriors on $256^3$8 from Licorice and Beorn versions of the same fiducial signal. The resulting biases, quantified as shifts in posterior means normalised by the larger posterior standard deviation, range from $256^3$9 to several $N = 1000$0 across five models, typically $N = 1000$1. The effect is strongest on $N = 1000$2: the Beorn-based posterior centres on a value more than 10% below truth despite near-identical mean temperatures, demonstrating that sub-noise agreement in aggregate quantities does not guarantee unbiased inference. The authors note that these model-error biases exceed those typically arising from switching inference methods while holding the simulator fixed. Because absolute parameter offsets remain modest ($N = 1000$3 dex in $N = 1000$4, $N = 1000$5 Gyr in $N = 1000$6), the SKA's sensitivity would resolve physically similar models but be vulnerable to this systematic scatter.

An important methodological caveat applies here: the study varies the inference target while keeping a single emulator trained on unmodified Licorice spectra, rather than training separate emulators per approximation — an expense beyond the paper's scope. The biases should therefore be read as order-of-magnitude estimates.

Limitations and open questions

The paper is candid about what remains unresolved. A portion of the power-spectrum difference is not explained by any tested approximation; the residual likely stems from the luminosity-binning preprocessing or effects not yet identified. Appendix results show that increasing the number of Beorn profiles from 250 to 20,000 changes the power spectrum by $N = 1000$715%, comparable to but smaller than the full Beorn–Licorice difference, and matching Licorice would require millions of profiles, eliminating Beorn's cost advantage. A single-source toy model confirms total deposited energy agrees closely but reveals profile-shape differences near the source, hypothesised to arise from Licorice's particles-to-cells mapping. Neither code can be declared more accurate: both embed uncertain modelling choices, and Licorice itself carries Monte Carlo noise and finite-resolution limitations. Generalisation to other codes (21cmFAST, Grizzly, C2-ray) is explicitly deferred to future work, though the demonstrated impact of the homogeneous-optical-depth approximation motivates a conservative expectation of significant disagreement with codes sharing it. Finally, whether additional nuisance parameters in fast codes, marginalised at inference time, could reconcile 1D and 3D posteriors while preserving computational savings remains an open question.

Conclusion

By processing identical Licorice source populations through Beorn, this work provides the first controlled comparison of X-ray heating between a 3D RT code and a 1D profile-painting code. Luminosities, mean temperatures, and global signals agree well, but temperature distributions, field topology, and 21-cm power spectra differ by $N = 1000$820–30%, sufficient to bias inferred astrophysical parameters by typically $N = 1000$9 at SKA-level noise. Individual common approximations — most notably the mean-density optical depth — are shown to produce discrepancies of this magnitude on their own. The ~20% power-spectrum agreement achieved here matches previous ionisation-focused comparisons, suggesting a robust floor for current methodology, and implies that interpreting future precision 21-cm data with heterogeneous simulation codes will require either substantially improved fast codes or explicit treatment of modelling error.