Thermal-bw speeds isotropic gamma-ray pair-production calculations

A new lightweight Python package calculates photon attenuation and the energy spectra of secondary electron–positron pairs in isotropic radiation fields, with its fitted blackbody approximation reaching a maximum relative error of 0.316% across its validated range. The software, called thermal-bw, targets a computationally demanding step in high-energy astrophysics: evaluating two-photon pair production, the process in which a gamma ray and a target photon create an electron and a positron. The package calculates total absorption, the differential energy distribution of either charged product and the pair source produced by an incident gamma-ray population. It assumes isotropic, unmagnetized and unpolarized radiation and is intended as a local interaction module for broader radiative-transfer, kinetic or Monte Carlo simulations rather than as a complete transport or cascade code.

The motivation comes from the structure of the calculation. In an isotropic field, the absorption coefficient requires integration of the Breit–Wheeler cross-section over both the energy of a target photon and the angle at which it meets a high-energy photon. Direct numerical evaluation is robust but can be slow, particularly when a model repeatedly evaluates many photon energies, temperatures or radiation-field components. For an ideal blackbody, the authors show that the opacity can be reduced to a universal function after its cubic temperature scaling is separated. Their new surrogate represents that function with six fitted parameters, uses only elementary operations and stores 48 bytes of float64 coefficients. A branchless construction makes the expression nonnegative for positive input energies and temperatures and avoids a need for piecewise conditionals in the hot path.

The dimensionless interaction variable in the approximation is defined as eta, equal to the product of the incident gamma-ray energy and the target thermal energy, each expressed relative to the electron rest mass. The surrogate was fitted over 0.01 <= eta <= 50 using 2,048 logarithmically spaced points, with fixed Gauss–Legendre quadrature employing 320 angular and 480 Planck-energy nodes as the reference. On 4,096 validation points selected between the fitting nodes, its median relative error was 0.225%, while its 99th-percentile and maximum errors were both 0.316%. A separate scrambled Sobol sample of 2,048 points produced the same maximum error. The largest logarithmic discrepancy was 0.00137 dex. Independent checks of the temperature reduction found that changing kT from 2.35 × 10⁻⁷ to 10³ keV altered the reduced coefficient by no more than 1.62 × 10⁻¹⁴ fractionally at fixed eta.

Outside the fitted interval, accuracy declines predictably. When the authors extended evaluation to 0.005 <= eta <= 100, the 99th-percentile and maximum relative errors rose to 5.30% and 5.48%, respectively, with the maximum at eta = 100. The package continues to permit such extrapolation, but its documentation identifies those results as outside the validated domain and recommends the general Gauss–Legendre route when high accuracy is required. The validated range corresponds to very different physical regimes depending on temperature: for a cosmic microwave background with kT near 2.35 × 10⁻⁷ keV, it spans incident gamma-ray energies from about 11.1 teV to 55.6 PeV; for a 1-keV blackbody, it runs from roughly 2.61 MeV to 13.1 GeV; and for a 30-keV blackbody, from about 87.0 keV to 435 MeV.

Speed is the approximation’s main practical advantage. In an archived benchmark performed on an AMD EPYC 9V74 processor, the surrogate evaluated up to roughly 10⁴ times faster than Gauss–Legendre quadrature and about 10⁷ times faster than adaptive quadrature for larger input arrays, with the speedup generally increasing as the arrays grew. A separate benchmark of the stored universal-function representations used 10⁶ eta values, four warm-up evaluations and 21 timed repetitions. The surrogate’s median evaluation time was 15.60 milliseconds, with a median absolute deviation of 0.18 milliseconds and a range from 15.22 to 17.15 milliseconds. Those measurements isolated evaluation of the function representation, excluding construction of eta and the common (kT)³ scaling factor. The reported timing used Python 3.13.5, NumPy 2.3.5 and SciPy 1.17.0 in a single process without worker threads or processes.

The package is more than a fast opacity estimator. Its isotropic routines can calculate the differential pair kernel and a normalized conditional distribution for the energy of one produced lepton. The electron and positron marginals are identical in the unpolarized isotropic treatment, and integration of the differential coefficient recovers the ordinary angle-averaged cross-section. The authors emphasize that assigning half of the incident photon’s energy to each particle is not generally equivalent to the calculated spectrum. For a 50-keV blackbody, the conditional distribution becomes increasingly unlike equal energy sharing as the incident gamma-ray energy rises. An example using an incident Eγ⁻² population interacting with that field produced an L1 difference of 0.154 between the normalized differential spectrum and an equal-energy-sharing approximation; among regions where the exact spectrum exceeded 1% of its peak, the largest fractional difference was 0.394.

Accuracy was checked through several routes. Four blackbody points agreed with 30-digit mpmath integration to a maximum relative difference of 5.7 × 10⁻⁹. For cutoff power-law and broken-power-law target spectra, adaptive and Gauss–Legendre opacity results agreed with independent 30-digit nested mpmath integrations to better than 4.1 × 10⁻⁹ in the tested cases. The differential kernel agreed with 50-digit calculations to better than 10⁻¹², including cases with substantially unequal photon energies, while the Böttcher–Schlickeiser analytical kernel and an independent centre-of-momentum-frame calculation recovered the same total cross-section to about 8 × 10⁻⁸ or better. Integrating the pair spectrum over several kinds of target field recovered the corresponding opacity, and the population routine conserved absorbed photon number and pair energy to about 2 × 10⁻⁶. These are numerical validation results for the software, not measurements of a physical process in an astrophysical environment.

The broader software supports blackbody, greybody, power-law, broken-power-law, tabulated, discrete-line, composite and user-callable target spectra. Known breaks, line energies and table knots are treated as integration boundaries. Callables with unknown narrow structure trigger increased quadrature order and comparison with adaptive integration before acceptance. Users can choose adaptive quadrature, fixed Gauss–Legendre integration or a cached method that interpolates the angle-averaged kernel while retaining integration over the supplied target spectrum. Differential calculations also include a centre-of-momentum-frame route to address cancellation in the analytical expression when the two photon energies differ greatly. The core requires Python 3.9 or later, NumPy 1.23 or later and SciPy 1.9 or later; optional Astropy support accepts physical quantities with units. Version 1.0.0 is distributed under an MIT license and uses Git for version control.

Comparisons with existing tools illuminate the package’s intended niche. Against a source-level reproduction of GAMERA’s ComputeAbsCoeff for a 50-keV blackbody, power law, broken power law and double-peaked spectrum, thermal-bw produced median absolute opacity differences of 2.23%, 2.09%, 2.36% and 1.92%, respectively, with maxima between 2.09% and 2.82%. When GAMERA’s spectrum interpolation and target-energy integration were retained but its approximate angle-averaged cross-section was replaced by thermal-bw’s numerical angular average, the median differences fell to 0.0127%, 0.0641%, 0.00739% and 0.0226%. The authors therefore attributed most of the 2–3% discrepancy to the approximate angular treatment rather than the spectral integration. A convention-adjusted comparison of the Breit–Wheeler cross-section itself with agnpy v0.5.1 found a maximum relative difference of 0.0 in float64 arithmetic across 2,000 logarithmically spaced test values. Direct end-to-end timing comparisons with AM3 and CRPropa were not possible because those packages perform broader calculations rather than the same isolated measurement.

The resulting library is designed to fit inside larger simulations that supply their own geometry, temporal evolution, particle cooling and escape. The same target object can be used to attenuate gamma rays and inject their secondary particles, while isotropic reference calculations can help test more elaborate angular or transport implementations. The package includes automated tests across Python 3.9–3.13 on Linux, macOS and Windows; a locally reproduced run passed 112 tests with 92.97% coverage. Its limits are equally explicit. It does not model source geometry, outgoing lepton angles, relativistic frame transformations, magnetic pair creation, cooling or cascade evolution, and it does not establish which numerical tolerance is appropriate for a particular scientific application. What it delivers is a compact, reproducible toolkit for local Breit–Wheeler opacity and pair spectra, with a sub-percent blackbody surrogate where its independently tested eta range applies and slower reference methods available elsewhere.

Subject of Research: Numerical modelling of Breit–Wheeler photon absorption and secondary pair spectra in isotropic radiation fields

Article Title: thermal-bw: A lightweight Python package for Breit–Wheeler opacity and secondary pair spectra in isotropic radiation fields

Article References: Harke, D. (2026). thermal-bw: A lightweight Python package for Breit–Wheeler opacity and secondary pair spectra in isotropic radiation fields. SoftwareX, 36, Article 103148. https://doi.org/10.1016/j.softx.2026.103148

Image Credits: AI Generated

DOI: 10.1016/j.softx.2026.103148

Keywords: thermal-bw, Breit–Wheeler process, gamma-ray opacity, pair production, blackbody radiation, Python software, high-energy astrophysics, radiative transfer, secondary particles, numerical validation, lightweight, Python