Numerical computation of isotropic Compton scattering

Compton scattering is involved in many astrophysical situations. It is well known and has been studied in detail for the past fifty years. Exact formulae for the different cross sections are often complex, and essentially asymptotic expressions have been used in the past. Numerical capabilities have now developed to a point where they enable the direct use of exact formulae in sophisticated codes that deal with all kinds of interactions in plasmas. Although the numerical computation of the Compton cross section is simple in principle, its practical evaluation is often prone to accuracy issues. These can be severe in some astrophysical situations but are often not addressed properly. In this paper we investigate numerical issues related to the computation of the Compton scattering contribution to the time evolution of interacting photon and particle populations. An exact form of the isotropic Compton cross section free of numerical cancellations is derived. Its accuracy is investigated and compared to other formulae. Then, several methods to solve the kinetic equations using this cross section are studied. The regimes where existing cross sections can be evaluated numerically are given. We find that the cross section derived here allows for accurate and fast numerical evaluation for any photon and electron energy. The most efficient way to solve the kinetic equations is a method combining a direct integration of the cross section over the photon and particle distributions and a Fokker-Planck approximation. Expressions describing this combination are given.

💡 Research Summary

The paper addresses a long‑standing problem in high‑energy astrophysics: the reliable numerical evaluation of the isotropic Compton scattering cross‑section for arbitrary photon and electron energies. While the exact Klein‑Nishina formula is well known, its direct implementation suffers from severe numerical cancellations and floating‑point overflow/underflow when the dimensionless energy ratio (r = \epsilon_\gamma/(\gamma m_ec^2)) spans many orders of magnitude. The authors first analyse the sources of these instabilities. In the traditional expression, logarithmic terms such as (\ln(1+2r)) are subtracted from nearly equal large numbers, and square‑root terms (\sqrt{1+2r}) become ill‑conditioned for both very small and very large (r). Consequently, standard double‑precision arithmetic can produce relative errors exceeding 10 % in the transition regime (10^{-4} \lesssim r \lesssim 10^{4}), which is precisely where many astrophysical plasmas operate.

To eliminate the cancellations, the authors reformulate the cross‑section in terms of two auxiliary functions, (A(r)) and (B(r)), that are constructed to remain bounded ((|A|,|B| \le 1)) across the entire range of (r). The new expression reads

\