Solving the transport equation by the use of 6D spectral methods in spherical geometry

We present a numerical method for handling the resolution of a general transport equation for radiative particles, aimed at physical problems with a general spherical geometry. Having in mind the computational time difficulties encountered in problems such as neutrino transport in astrophysical supernovae, we present a scheme based on full spectral methods in 6d spherical coordinates. This approach, known to be suited when the characteristic length of the dynamics is much smaller than the domain size, has the potential advantage of a global speedup with respect to usual finite difference schemes. An analysis of the properties of the Liouville operator expressed in our coordinates is necessary in order to handle correctly the numerical behaviour of the solution. This reflects on a specific (spherical) geometry of the computational domain. The numerical tests, performed under several different regimes for the equation, prove the robustness of the scheme: their performances also point out to the suitability of such an approach to large scale computations involving transport physics for mass less radiative particles.

💡 Research Summary

The paper presents a novel numerical framework for solving the six‑dimensional transport equation that governs the evolution of massless radiative particles such as photons and neutrinos. Recognizing that traditional finite‑difference (FD) methods become prohibitively expensive in high dimensions, the authors adopt a global spectral approach formulated in spherical coordinates for both physical space (r, θ, φ) and momentum space (ε, θₚ, φₚ).

First, the Liouville operator is rewritten in these coordinates, exposing singular behavior at the origin (r = 0). To avoid numerical instability, the authors propose a hybrid discretization: Chebyshev expansions are used for angular and energy variables, while the radial direction can be treated either with a Chebyshev basis or with a conventional second‑order finite‑difference stencil. This flexibility allows the method to handle sharp radial gradients without sacrificing spectral accuracy elsewhere.

The collision term is modeled in two physically relevant limits. In the coherent‑scattering regime, energy is conserved (δ(ε‑ε′)) and the differential cross‑section depends only on the cosine of the scattering angle. For photons this reduces to the Thomson kernel; for neutrinos it takes a simple weak‑interaction form. In the Fokker‑Planck approximation, valid for low‑energy particles in a hot plasma, the integral scattering operator is replaced by drift‑diffusion operators in energy space, introducing two characteristic timescales: an isotropisation time τ_is and a much longer “bosonisation” time τ_B. The authors derive the corresponding reduced equations, including a Kompaneets‑type equation for the photon distribution, and demonstrate analytically that particle number is conserved.

Implementation details are given for a highly parallel code. Six‑dimensional tensors are distributed across MPI processes, while OpenMP threads accelerate local matrix‑vector products. Time integration employs operator splitting: the free‑streaming Liouville term is advanced with an explicit spectral propagator, and the collision term is treated implicitly or semi‑implicitly to maintain stability.

A suite of benchmark problems validates the algorithm. (1) Free propagation of a uniform distribution confirms that the spectral scheme preserves the analytic solution to machine precision. (2) Pure coherent scattering tests the angular redistribution and particle‑number conservation. (3) The Fokker‑Planck test reproduces the expected rapid isotropisation of the angular distribution and the slower evolution toward a Bose‑Einstein spectrum, matching analytical predictions for τ_is ≪ τ_B. (4) An application to a rotating neutron star combines a diffusion approximation in the dense core with an exact transport solution in the outer region, using the telegraph equation to bridge the two regimes. The telegraph formulation, which limits signal speed to c/√3, is shown to be more appropriate than pure diffusion when matter velocities approach relativistic values.

Performance analysis shows that, for a target error of 10⁻⁶, the spectral method requires roughly five times fewer grid points per dimension than a second‑order FD scheme, leading to an overall reduction in degrees of freedom by a factor of 5⁶ ≈ 1.6 × 10⁴. Consequently, the total floating‑point operation count drops by about six orders of magnitude. On a single 2.5 GHz core the full six‑dimensional test runs complete in ~30 minutes, and the code exhibits near‑linear scaling up to thousands of cores on modern HPC clusters.

The authors also provide three appendices: (A) an explicit particle‑conserving formulation of the transport equation; (B) a fully spectral solution for the spherically symmetric (1‑D) case, demonstrating that high‑frequency modes are accurately captured; and (C) derivations of the diffusion and telegraph equations, arguing that the latter is preferable for relativistic flows because it respects causality while remaining straightforward to implement within the spectral framework.

In summary, the work demonstrates that a six‑dimensional spectral method in spherical coordinates can achieve high accuracy, robust conservation properties, and dramatic computational savings compared with conventional finite‑difference approaches. This makes it a compelling tool for large‑scale astrophysical simulations involving neutrino or photon transport, such as core‑collapse supernovae, neutron‑star mergers, and high‑energy X‑ray sources, where both complex geometry and relativistic matter motions are essential.