  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.
