EigenFlux: a stable multi-stream radiative transfer method for strongly scattering media

Radiative transfer in strongly scattering media remains computationally challenging, particularly for systems with highly asymmetric phase functions and large optical depths, where conventional discrete-ordinate approaches may suffer from numerical instability, slow convergence, or loss of accuracy. We present EigenFlux, a multistream radiative transfer framework based on eigenmode decomposition, natural-reflectance stabilization, and flexible mesh-based angular discretization. Unlike previous approaches relying primarily on global polynomial expansions, EigenFlux permits localized resolution of phase functions with strongly forward scattering or significant backward-scattering components while preserving flux conservation and numerical stability. The method is evaluated across 957 test cases spanning asymmetry factors from g=-0.998867 to g =0.998867 and single-scattering albedos from ω =0.0014660 to ω =0.9999995, including extreme multiple-scattering regimes that are difficult for conventional solvers. Comparisons with DISORT show 748 times better accuracy and 10 times faster execution for number of streams greater or equal to 168. EigenFlux maintained stable solutions for asymmetry factors exceeding ∣g∣>0.99, although stability is weaker for absorption levels less than 0.10. Analysis of the eigenspectrum reveals the emergence of asymptotic diffuse transport regimes in optically thick systems (thick/deep enough for the reflectance to reach its natural limit) and persistent direct-beam structure in semi-transparent media. These results suggest that EigenFlux provides a stable and flexible framework for radiative transfer calculations for atmospheres, snow and ice, ocean optics, pigments and coatings, remote sensing, and graphics rendering, particularly in cases involving extreme scattering asymmetry.