Observations and magnetohydrodynamic simulations of solar and stellar atmospheres reveal an intermittent behavior or steep gradients in physical parameters, such as magnetic field, temperature, and bulk velocities. The numerical solution of the stationary radiative transfer equation is particularly challenging in such situations, because standard numerical methods may perform very inefficiently in the absence of local smoothness. However, a rigorous investigation of the numerical treatment of the radiative transfer equation in discontinuous media is still lacking. The aim of this work is to expose the limitations of standard convergence analyses for this problem and to identify the relevant issues. Moreover, specific numerical tests are performed. These show that discontinuities in the atmospheric physical parameters effectively induce first-order discontinuities in the radiative transfer equation, reducing the accuracy of the solution and thwarting high-order convergence. In addition, a survey of the existing numerical schemes for discontinuous ordinary differential systems and interpolation techniques for discontinuous discrete data is given, evaluating their applicability to the radiative transfer problem.

The numerical computation of reliable and accurate Stokes profiles is of great relevance in solar physics. In the synthesis process, many actors play a relevant role: among them the formal solver, the discrete atmospheric model, and the spectral line. This paper tests the performances of different numerical schemes in the synthesis of polarized spectra for different spectral lines and atmospheric models. The hierarchy between formal solvers is enforced, stressing the peculiarities of high-order and low-order formal solvers. The density of grid points necessary for reaching a given accuracy requirement is quantitatively described for specific situations.

Efficient numerical approximation of the polarized radiative transfer equation is challenging because this system of ordinary differential equations exhibits stiff behavior, which potentially results in numerical instability. This negatively impacts the accuracy of formal solvers, and small step-sizes are often necessary to retrieve physical solutions. This work presents stability analyses of formal solvers for the radiative transfer equation of polarized light, identifies instability issues, and suggests practical remedies. In particular, the assumptions and the limitations of the stability analysis of Runge-Kutta methods play a crucial role. On this basis, a suitable and pragmatic formal solver is outlined and tested. An insightful comparison to the scalar radiative transfer equation is also presented.

The discussion regarding the numerical integration of the polarized radiative transfer equation is still open and the comparison between the different numerical schemes proposed by different authors in the past is not fully clear. Aiming at facilitating the comprehension of the advantages and drawbacks of the different formal solvers, this work presents a reference paradigm for their characterization based on the concepts of order of accuracy, stability, and computational cost. Special attention is paid to understand the numerical methods belonging to the Diagonal Element Lambda Operator family, in an attempt to highlight their specificities.

When integrating the radiative transfer equation for polarized light, the necessity of high-order numerical methods is well known. In fact, well-performing high-order formal solvers enable higher accuracy and the use of coarser spatial grids. Aiming to provide a clear comparison between formal solvers, this work presents different high-order numerical schemes and applies the systematic analysis proposed by Janett et al. (2017), emphasizing their advantages and drawbacks in terms of order of accuracy, stability, and computational cost.