Discrete Ordinates Approximation

• Discrete ordinates approximation is a method that discretizes the angular variable using fixed directions and quadrature weights, converting continuous integrals into tractable algebraic or differential equations.
• It utilizes various spatial discretization techniques and stabilization methods—such as DODSD and weak Galerkin—to enhance accuracy, control numerical oscillations, and mitigate ray effects.
• Ongoing research focuses on hybrid schemes, low-rank tensor approaches, and Monte Carlo integrations to accelerate convergence and extend DOM’s applicability to complex geometries and multi-physics problems.

The discrete ordinates approximation, commonly referred to as the DOM or Sₙ method, is a foundational approach for numerically solving the radiative transfer equation (RTE), linear Boltzmann equation, and related linear kinetic equations. By discretizing the angular variable—replacing integrals over continuous directions with sums over a finite set of discrete directions with associated quadrature weights—DOM reduces the original integro-differential system to a large but tractable system of coupled differential or algebraic equations. The method forms the backbone of deterministic radiation transport calculations in domains as varied as nuclear engineering, atmospheric science, astrophysics, and biomedical imaging. Its flexibility has led to continual development of improved spatial and angular discretizations, acceleration schemes, and hybridizations with Monte Carlo algorithms.

1. Angular Discretization and the Discrete Ordinates Principle

The discrete ordinates method transforms an angular integral such as

$\int_{S^{d-1}} \psi(\mathbf{x}, \Omega)\, d\Omega$

into a quadrature sum

$\sum_{n=1}^{N} w_n\, \psi(\mathbf{x}, \Omega_n)$

where each direction $\Omega_n$ is one of a finite set of prescribed unit vectors and $w_n$ is the associated quadrature weight. For example, in slab geometry, the integration over $\mu\in [-1,1]$ is replaced by $N$ Gauss-Legendre points. In full three-dimensional geometries, Level Symmetric, Gauss-Legendre, or product quadrature sets are commonly used (Ganapol, 2014, Singh, 2022, Ma et al., 20 Mar 2025).

After angular discretization, the RTE becomes a system of $N$ transport equations, each for the intensity along a discrete direction: $$\Omega_n \cdot \nabla \psi_n(\mathbf{x}) + \sigma_t(\mathbf{x}) \psi_n(\mathbf{x}) = \sum_{m=1}^N w_m \sigma_s(\mathbf{x}) \Phi(\Omega_m, \Omega_n) \psi_m(\mathbf{x}) + Q_n(\mathbf{x})$$ with boundary or initial conditions specified on inflow boundaries (Wang et al., 2016, Singh, 2022).

This approximation's accuracy and computational cost are governed by the quadrature order $N$. Low $N$ leads to the so-called ray effect—a major artifact discussed below.

2. Numerical Methods, Stability, and Error Analysis

After angular discretization, the spatial part of the resulting system is typically discretized using finite difference, finite volume, or finite element techniques.

A class of stabilized finite element approaches, such as the discrete-ordinate discontinuous-streamline diffusion (DODSD) (Wang et al., 2016) and weak Galerkin (DOWG) (Singh, 2022) methods, is designed to address the hyperbolic nature of each directionally discretized equation. The DODSD method introduces an artificial streamline diffusion term, replacing standard test functions $v$ with $v + \delta\, (\omega\cdot\nabla v)$, which helps control numerical oscillations and improves error estimates, especially in the directional derivative norms. A crucial result for these methods is a stability bound such as: $$||| \{u_h^l \} |||^2 \leq 3 a_h( \{u_h^l \}, \{u_h^l \} )$$ where $a_h$ is the discrete bilinear form and $||| \cdot |||$ is an appropriately defined norm incorporating both $L^2$ and gradient-like terms (Wang et al., 2016).

Error bounds for these methods typically have the form

$$||| \{u^l\} - \{u_h^l\} ||| \leq C_1 h^{\min\{r,k\} + 1/2} \left( \sum_l \|u^l\|^2_{H^{r+1}(X)} \right)^{1/2}$$

where $h$ is the spatial mesh size, $k$ is the polynomial degree, and $r$ the regularity of the true solution.

The weak Galerkin method achieves similar results without requiring tuning of penalty parameters, offering robustness particularly on general meshes and complex domains (Singh, 2022).

3. Acceleration, Ray Effect Mitigation, and Hybrid Schemes

Ray effects are a direct result of approximating the angular integral by a fixed sparse set of directions: the solution exhibits spurious oscillations along the rays defined by the discrete ordinates. Several strategies have been developed to mitigate this problem:

(a) Quasi-Random and Random Ordered Methods

The quasi-random discrete ordinates (QRDOM) (Konzen et al., 2023, Konzen et al., 2023) and random ordinate method (ROM) (Li et al., 17 Jul 2024) replace the fixed quadrature set with sampling-based angular discretizations. QRDOM uses low-discrepancy sequences (e.g., reverse Halton sequence) to sample directions, smearing the regularity that leads to ray effects, while ROM employs random sampling of direction in angular subcells.

Key formulas for the ROM error and bias decay are

$$\mathbb{E}\|\phi^{(\xi)} - \phi\| \leq C n^{-3/2} (\log n)^{1/2}$$

$$\|\mathbb{E} \phi^{(\xi)} - \phi\| \leq C \lambda n^{-3} \log n$$

where $n$ is the number of angular cells (Li et al., 17 Jul 2024). This demonstrates that the ensemble-averaged solution has significantly improved convergence and greatly reduces the ray effect while maintaining computational cost.

(b) Implicit and Low-Rank Methods

The AREPO-IDORT scheme (Ma et al., 20 Mar 2025) implements an implicit discrete ordinates transport method on an unstructured moving mesh, fully coupled to an MHD solver. The implicit solver enables large timesteps (not limited by cell light-crossing time), while careful design of angular quadratures preserves directional fidelity and reduces artifacts in both optically thin and thick domains. By evolving all angles simultaneously with local time-stepping and interfacing to a Riemann solver, AREPO-IDORT enables accurate multi-scale radiation-MHD simulations.

Low-rank tensor approximations (Peng et al., 2022) further reduce the computational complexity of time-dependent multi-dimensional Sₙ solvers by approximating the flux as a low-rank decomposition in space and angle, preserving the upwind sweeping structure required for the transport solver. For appropriate ranks $r$, these methods demonstrate that memory usage and runtime can be reduced by over one order of magnitude compared to full-rank approaches, with minimal accuracy loss.

(c) Hybrid MC–Discrete Ordinates Schemes

Hybridization of DOM with Monte Carlo, as in the MC–DG hybrid (Krotz et al., 2023), exploits the strengths of each method: Monte Carlo efficiently simulates the highly directional "uncollided" flux (which is where DOM ray effects are most pronounced), while DOM handles the smoother "collided" component. The resulting approach achieves higher accuracy and substantially reduced computational complexity compared to pure DOM in challenging problems, especially for sharp wavefronts and deep penetration geometries.

4. Analytical Solution Strategies and Convergence Acceleration

In one-dimensional homogeneous and layered media, analytical and semi-analytical approaches improve both stability and precision.

The response matrix (RM/DOM) formulation (Ganapol, 2014) recasts the solution in terms of matrix hyperbolic functions: $$y(t) = H(t)\, y(T_0) + [H(T_0 - t) - H(t)]\, (\text{source term}), \quad H(t) = \sinh(V_A t)$$ where $V_A$ is assembled from the eigenvalues/vectors of the angular operator. This avoids the catastrophic cancellation associated with exponentials of widely separated eigenvalues, making the method robust for highly conservative or near-conservative cases.

Convergence of integrated quantities (e.g., reflectance, transmittance) as a function of quadrature order is accelerated using Wynn-epsilon ($\epsilon$) algorithms: $$\epsilon_{-1}(N) = 0, \quad \epsilon_0(N) = R_f(N), \quad \epsilon_{k+1}(N) = \epsilon_{k-1}(N+1) + \frac{1}{\epsilon_k(N+1) - \epsilon_k(N)}$$ This enables computation of slab quantities to $10^{-12}$ relative error with modest $N$, even for thin slabs and challenging phase functions.

The same analytical approach can be extended, via the "star product" formulation of response matrices, to heterogeneous slab stacks with different optical properties—enabling recursive, stable solution assembly for complex layered systems (Ganapol, 2014).

5. Extensions to Heterogeneous Media, Complex Geometries, and Advanced Physics

DOM is highly extensible. Layered and heterogeneous systems are efficiently handled via response matrix or star product constructions (Ganapol, 2014, Leonardo et al., 2021), with adjoint and importance solutions computed via specialized nodal methods (e.g., spectral nodal adjoint response matrix).

For complex 2D and 3D geometries, meshless methods combining DOM angular discretization with moving least squares (MLS) spatial approximation have demonstrated robust accuracy (up to $2\%$ relative error in challenging 3D enclosures) without resorting to mesh generation—a practical advantage for radiative transfer in complex domains (Wang et al., 2018).

Novel DOM formulations, such as the hybrid discrete H$^T_N$ scheme (McClarren et al., 2021), combine the P$_N$ (spherical harmonics) and S$_N$ (standard DOM) approaches by representing the angular variable as a piecewise Legendre expansion within angular bands. This enables better control of oscillations and negative fluxes typical in P$_N$, while mitigating the ray effects of S$_N$. The H$^T_N$ model produces a hyperbolic system for all $T \geq 1$, $N \geq 0$ and preserves asymptotic diffusion limits.

For nonclassical transport, spectral expansion in Laguerre polynomials (in the free-path variable) combined with Sn methods enables efficient solution of transport equations where memory effects are important, though with limitations related to convergence of the Laguerre moments (requiring $\sigma_t > (\sqrt{3}/6)$) and numerical stability at high truncations (Moraes et al., 2021).

6. Practical Implementations and Benchmarks

Large-scale practical systems, including X-ray tomography, atmospheric remote sensing, and reactor analysis, employ advanced DOM-based codes. For example:

• The DOCTORS code (Norris et al., 2018, Norris et al., 2019) implements DOM discretization in space, energy, and angle, with first-collision source treatment and GPU acceleration, enabling near real-time CT dose calculations. RMSD with respect to Monte Carlo is ~2.87% when protocol parameters (mesh, angular and Legendre order) are chosen appropriately.
• Remote sensing and optical tomography applications, such as those using the SHDOM framework (Levis et al., 2015), combine discrete ordinates with spherical harmonics expansions and efficient optimization-based inversion, enabling 3D recovery of extinction fields in clouds, validated on $>10^5$ grid points with errors ($\lesssim 1\%$ mean).
• Benchmark comparisons of numerical schemes (LD, step-characteristic, exponential-discontinuous, etc.) reveal the trade-offs between accuracy, positivity, and computational cost in practical slab transport (Roberts, 14 Feb 2024).

Testing against exact or semi-analytic solutions (e.g., for Henyey–Greenstein scattering, cloud fields with spatial heterogeneity) demonstrates DOM's capacity for high-precision (7–8 digit) benchmarks with proper formulation and acceleration techniques (Ganapol, 2014).

7. Limitations, Research Trends, and Future Directions

Despite its widespread use and continual innovation, the discrete ordinates approximation exhibits several persistent challenges:

• Ray effects remain a key limitation for standard DOM in low-quadrature or non-smooth solution regimes, addressed through stochastic angular sampling, hybrid MC–DOM schemes, and domain-adapted quadrature generation (Konzen et al., 2023, Li et al., 17 Jul 2024, Krotz et al., 2023).
• Acceleration for strongly forward-peaked (nearly singular) scattering kernels requires advanced preconditioning, e.g., Fokker–Planck synthetic acceleration, to efficiently damp higher-order angular errors (Patel et al., 2019).
• The curse of dimensionality continues to drive the exploration of low-rank tensor approaches and meshless spatial solvers (Peng et al., 2022, Wang et al., 2018).
• Adaptive and hybrid angular–spatial discretizations that automatically balance accuracy and computational cost are emerging as viable strategies, especially for multi-physics simulations (e.g., radiation–MHD) (Ma et al., 20 Mar 2025).

A plausible implication is that, with the continued development of implicit solvers, randomized angular quadrature, and integration into massively parallel architectures, DOM will remain core to exascale simulation and high-resolution, multi-scale modeling in physical and engineering sciences. The evolving landscape of hybrid schemes, data-driven acceleration, and meshless models also signals increased generality and robustness for future DOM-based applications.