Discontinuities in numerical radiative transfer
^{1}
Istituto Ricerche Solari Locarno (IRSOL), 6605 LocarnoMonti, Switzerland
email: gioele.janett@irsol.ch
^{2}
Seminar for Applied Mathematics (SAM) ETHZ, 8093 Zurich, Switzerland
Received:
30
July
2018
Accepted:
27
November
2018
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 firstorder discontinuities in the radiative transfer equation, reducing the accuracy of the solution and thwarting highorder 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.
Key words: radiative transfer / methods: numerical / polarization
© ESO 2019
1. Introduction
The presence of sharp variations in the physical parameters of solar and stellar atmospheres is anything but irrelevant in the formation of spectroscopic and spectropolarimetric signals. For instance, it is well known that the propagation of shock waves can generate a large variety of important features in observed spectra (Mihalas 1981). Moreover, the presence of transition layers or boundaries associated with a jump in magnetic field, in temperature, and in bulk velocities significantly modifies spectral line profiles.
In many astrophysical applications, one must routinely solve the radiative transfer equation. Very common examples are: iterative solutions of nonlinear radiative transfer problems in nonlocal thermodynamic equilibrium (nonLTE) conditions, inversions of spectral line profiles, and massive threedimensional (3D) radiative magnetohydrodynamic simulations (RMHD) of stellar atmospheres.
The radiative transfer equation seldom has solutions that can be expressed in an analytical form, and therefore it is common to seek approximate solutions in terms of numerical schemes. However, standard numerical methods rely on smoothness assumptions regarding the functions to be integrated and these methods may become inaccurate or inefficient in the presence of discontinuities. In particular, the crossing of a discontinuity might introduce important local errors in the numerical solution, which negatively affect the accuracy of the global problem.
The relevance of discontinuities has already been recognized in the radiative transfer community. For example, standard nonLTE radiative transfer codes barely converge when applied to data from 3D RMHD simulations that contain discontinuities and steep gradients (Steiner et al. 2016). Some efforts have already been exercised to handle such discontinuities: Tscharnuter (1977) proposed a numerical method for computing the radiative transfer across a shock front with spherical symmetry. Afterwards, Mihalas (1981) designed a numerical scheme for solving the equation of transfer in discontinuous media. Auer & Paletou (1994) remarked the necessity to avoid negative overshoots when interpolating physical quantities. For this reason, Auer (2003) suggested the use of Hermite interpolants and Bézier curves in the formal solution and, subsequently, Ibgui et al. (2013) made use of monotonic cubic Hermite polynomials.
Some attention to discontinuities has also been paid in the context of the radiative transfer of polarized light. Steiner (2000) and Müller et al. (2002) studied the role played by discontinuities in fluid velocities and in the magnetic field in the formation of asymmetric Stokes profiles. The SIRJUMP inversion code by Louis et al. (2009) is able to infer possible discontinuities in the physical quantities along the line of sight (del Toro Iniesta & Ruiz Cobo 2016). More recently, de la Cruz Rodríguez & Piskunov (2013) and Štěpán & Trujillo Bueno (2013) applied Bézier interpolants to control abrupt changes in the absorption and emission quantities, whereas Steiner et al. (2016) proposed a numerical scheme based on monotonic reconstruction that allows for discontinuities at the boundary of each numerical cell of the atmospheric model.
However, a rigorous mathematical and numerical investigation of the problems arising in the integration of the discontinuous radiative transfer equation is still lacking and deeper understanding of this topic is certainly required. This work clearly states the issues introduced by discontinuities in the context of the radiative transfer equation. The paper is organized as follows: Sects. 2 and 3 introduce differential systems with a discontinuous righthand side and interpolations of discontinuous functions, respectively, explaining the limitations of standard convergence analyses for these kinds of problems. Section 4 explores the role of discontinuities in the radiative transfer equation. Section 5 exposes specific numerical tests that highlight the inefficient performances of standard methods when dealing with discontinuities in the radiative transfer equation. Sections 6 and 7 briefly summarize the existing numerical methods for the treatment of discontinuous ordinary differential equations (ODEs) and the interpolation techniques for discontinuous discrete data. Particular attention is paid to their suitability to the radiative transfer equation in discontinuous media. Finally, Sect. 8 provides remarks and conclusions, with a view on future work.
2. Ordinary differential equations with a discontinuous righthand side
Discontinuous differential systems frequently appear in various fields, such as chemical kinetics (e.g., Landry et al. 2009), biological systems (e.g., Gouzé & Sari 2002), and various engineering disciplines (e.g., Malmborg & Bernhardsson 1997; Acary & Brogliato 2008). Moreover, most of the studies on discontinuous ODEs rely on numerical methods and a broad literature attests to the importance of this topic: from the pioneering works by Chartres & Stepleman (1972) and Mannshardt (1978) to the surveys by Calvo et al. (2008) and Dieci & Lopez (2012).
2.1. Impact of a discontinuity on the numerical solution
The numerical integration of discontinuous differential systems is challenging, because standard numerical schemes may perform very inefficiently in the absence of local smoothness. It is known that the crossing of a discontinuity introduces significant local errors, which negatively affect global errors and might even yield a convergence order breakdown.
Consider a numerical method of order p applied to the scalar initial value problem (IVP)
Standard convergence analysis would conclude that the numerical scheme yields a global error that scales as 𝒪(h^{p}), where h denotes the step size. Such an analysis relies on the assumption that the righthand side f of the IVP (1) (and consequently y) is sufficiently smooth. More precisely, the function f must have p + 1 continuous derivatives. However, a different local truncation error analysis is mandatory if discontinuities occur in f or its derivatives. By definition, a discontinuity in Eq. (1) is of order q ≥ 1 if f contains a finite jump in at least one of its partial derivatives of order q − 1 and it has continuous derivatives through order q − 2 (Gear & Osterby 1984). This produces a qthorder discontinuity in y at the discontinuity location.
Numerical analysis shows that if q ≥ p + 1 the numerical method will remain of order p and one may not even notice the discontinuity. By contrast, if q ≤ p, the discontinuity introduces the dominant term 𝒪(h^{q}) in the local error, which can reduce the global error order. For instance, Mannshardt (1978) observed that a RungeKutta method remains convergent after having crossed a firstorder discontinuity, but only with order 1.
In many practical problems, the righthand side f of the IVP (1) is known at a discrete set of points only. In such cases, the numerical grid must be considered fixed a priori and there is not a well defined concept of discontinuity in f. It is difficult to tell whether the underlying model that gave rise to the discrete data is continuous, or rather contains jumps. In a discrete world, all changes are jumps by definition and data are purely discontinuous even if the underlying model is not. Therefore, there is no concrete difference between a discontinuity and a steep gradient for a single set of discrete data that is not able to resolve the local feature as shown in Fig. 1. The difference is only observed if one considers a sequence of discrete grids with a finer and finer sampling. Computing the maximum of the differences between adjacent data points along the sequence of increasingly refined grids, one recognizes two possible cases: either the underlying function is continuous and the value approaches zero, or it contains a discontinuity and the value approaches the size of the jump.
Fig. 1. The hyperbolic tangent (top panel, continuous blue) and the modified Heaviside (bottom panel, continuous blue) functions are sampled in the interval x ∈ [−1, 1] with 10 (red circles) and 40 (green dots) grid points. No difference between the coarser (red) samplings is visible, while this is not true for the finer (green) samplings. 

Open with DEXTER 
Practical radiative transfer problems often involve rapidly varying functions which actually do not contain any mathematical discontinuities. Section 5.1 numerically explores such issues by considering both a discontinuous atmospheric model and one containing a steep gradient.
2.2. Example: trapezoidal method order reduction
This section gives a concrete example of order reduction for the (implicit) trapezoidal method, but the following local error analysis can be adapted to other numerical methods. For instance, Mannshardt (1978) presented the case of a general onestep method, whereas Gear & Osterby (1984) analyzed predictorcorrector methods.
For simplicity, consider the IVP
where f : ℝ → ℝ and the discretization of the time interval [0, T] given by {t_{k}} with k = 0, …, N. Let ξ ∈ [t_{k}, t_{k + 1}] and assume both f_{1} and f_{2} at least threetimes continuously differentiable. For f_{1}(ξ)≠f_{2}(ξ), the IVP (2) shows a firstorder discontinuity. The trapezoidal method applied to the IVP (2) in the interval [t_{k}, t_{k + 1}] reads
where y_{k} and y_{k + 1} are the numerical approximations of the exact values y(t_{k}) and y(t_{k + 1}), respectively. The local truncation error is defined by
while the global error, defined as
represents the accumulation of local errors over all the steps. The Taylor expansion of the local exact solution gives
Making use of the additional Taylor expansions
one recovers by direct calculations the expression
where h = t_{k + 1} − t_{k}. This leads to the formula
where
are the first and secondorder jumps, respectively, at the discontinuity point ξ. The first and second terms in Eq. (4) contribute to the global error with 𝒪(h) and 𝒪(h^{2}) terms, respectively. These contributions arise from the jump and cannot be improved by using some higherorder method. After having crossed a firstorder discontinuity, that is, f_{1}(ξ)≠f_{2}(ξ), the trapezoidal method is therefore only firstorder convergent. Here, one assumes a finite number of jumps, otherwise the 𝒪(h) contributions will accumulate, thwarting convergence. Discontinuities of second order or greater (which imply f_{1}(ξ) = f_{2}(ξ)) do not affect the order of convergence of the trapezoidal method.
3. Interpolations of discontinuous functions
Consider a function y(x) and a discrete set of points x_{0} < x_{1} < … < x_{n} and assume that y(x_{0}),…,y(x_{n}) are given. The interpolation problem is to find a function p(x) that satisfies interpolation requirements
The simplest way is to connect the given points with straight lines. However, when seeking for higher accuracy, it is usual to look for a function p(x) that is a polynomial or a piecewisesmooth polynomial of higher degree. Alternatively, the function p(x) can be a trigonometric function, a rational polynomial function, and so on.
3.1. Polynomial interpolation
In the problem of polynomial interpolation, one seeks a polynomial p(x) that satisfies the interpolation condition (5). Limiting the degree of p(x) to be ≤n, one obtains precisely one single interpolant p_{n}(x) that satisfies the interpolation conditions. A very common strategy is to construct this interpolant in terms of Lagrange polynomials, namely
where the Lagrange basis polynomials, ℓ_{i}(x), given by
satisfy the relation ℓ_{i}(x_{j}) = δ_{ij}, with δ_{ij} being the Kronecker delta
The interpolant p_{n}(x) and the interpolated function y(x) in Eq. (6) satisfy the condition (5), that is, they agree with each other at the interpolation points. In general, there is no reason to expect them to be close to each other elsewhere. Nevertheless, one can estimate the difference between them, the socalled interpolation error.
Let Π_{n} denote the space of polynomials of degree ≤n, and let C^{n + 1}[a, b] denote the space of functions that have n + 1 continuous derivatives on the interval [a, b]. Suppose y(x)∈C^{n + 1}[a, b] and let p_{n}(x)∈Π_{n} be the unique polynomial that interpolates y(x) at the n + 1 distinct points x_{0} < x_{1} < … < x_{n} ∈ [a, b]. Then, for all x ∈ [a, b] there is a ξ ∈ (a, b), such that
Defining h = max_{j} x_{j + 1} − x_{j}, one obtains
This means that, by assuming y(x) smooth enough, the error decreases to zero as 𝒪(h^{n+1}) and the interpolation has an order of accuracy n + 1. However, the smoothness assumption of y(x) is not always fulfilled and polynomial interpolations may become inaccurate in the presence of discontinuities.
Moreover, the error of the interpolation does not necessarily decrease by increasing the order of the polynomial. In fact, even though h^{n + 1} may go to zero for n → ∞, the term y^{(n + 1)} can grow rapidly, preventing convergence in Eq. (8). Equidistant interpolation of the Runge’s function (the socalled Runge phenomenon) is a striking example of this.
3.2. Impact of discontinuities
Error estimates (7) and (8) assume some smoothness of y(x) and the order of accuracy of interpolations decreases for less regular functions. Suppose that y(x)∈C^{s}[a, b] with s < n. Then, Eq. (7) must be replaced by
where is any subset of . By the same arguments that led up to Eq. (8), one obtains
that is, the order of accuracy of the polynomial interpolation is reduced to s.
Moreover, it is known that standard secondorder or higher interpolations are oscillatory near discontinuities and such oscillations might lead to numerical inaccuracy and even to numerical instability in nonlinear problems (e.g., Richards 1991; Shu 1998). This implies, first, an over/undershoot in p_{n}(x) (n > 1) whenever the function y(x) has a jump discontinuity and, second, the failure of higherorder approximations to remove the overshoot. This behavior (similar to the Gibbs phenomenon in Fourier series) is known as Gibbs phenomenon. For example, Zhang & Martin (1997) showed that cubic spline interpolations on uniform meshes always oscillate near a jump discontinuity. For these reasons, it makes little sense (and could even be noxious) to use standard highorder interpolation (n > s) when the smoothness of y(x) is not guaranteed.
3.3. Interpolatory quadratures
Consider a function y(x) and assume that a discrete set of points x_{0} < x_{1} < … < x_{n} and the data values y(x_{0}),…,y(x_{n}) are given. A systematic approximation of the integration of y(x) in an interval [ã, b̃] is commonly called quadrature and its formula is usually given in the form
where ω_{i} are the quadrature weights. A usual way to derive a quadrature formula is through the socalled interpolatory quadrature: one makes an assumption for the polynomial form of the integrand in the integration interval and analytically solves the integral. For instance, the Lagrange interpolation given by Eq. (6) yields
where the integrals of the Lagrange basis polynomials provide the quadrature weights. The quadrature error, i.e., the difference between the lefthand and the righthand sides of Eq. (10), decreases to zero as 𝒪(h^{n+2}). Moreover, the quadrature is exact for all polynomials of degree ≤n, that is, for all y(x)∈Π_{n}.
However, Sect. 3.2 shows that a polynomial of degree n guarantees 𝒪(h^{n+1}) accuracy only if y(x)∈C^{n + 1} in the integration interval. Consequently, interpolatory quadratures suffer from the same orderreduction issue of interpolations when facing discontinuities. Moreover, the presence of overshoots in the interpolated integrand reduces the accuracy of the quadrature and could even yield negative results in the quadrature of positive integrands. For this reason, the use of standard highorder interpolatory quadratures is not encouraged when the smoothness of y(x) is not guaranteed.
4. Discontinuities in radiative transfer
Standard radiative transfer calculations usually assume smooth variations in the radiation field and in the atmospheric physical parameters. The incidence of discontinuities is often neglected, without considering whether this assumption is justified or not. This section inquires into the role of discontinuities in the radiative transfer equation, focusing on the possible numerical issues.
4.1. Discontinuities in solar and stellar atmospheres
The plasma of solar and stellar atmospheres can be highly dynamic, with its properties fluctuating sharply over small scales. For instance, jumps in the magnetic field frequently appear: a typical example is the magnetopause at the boundary of a flux tube, where the magnetic field of the tube is separated from the surrounding fieldfree plasma by a thin transition layer. Another example is given by magnetic canopies, that is, horizontally extending magnetic fields overlying fieldfree regions of internetwork cells and separated from it by a thin transition layer (Steiner 2000). Abrupt variations in temperature can also be present: for example, the jump in temperature of three orders of magnitude that takes place in the solar transition region. Clouds of cool material embedded in the chromosphere or corona are usually associated with jumps in bulk velocities as well. For instance, a magnetic flux tube interacting with convective motions shows a downwardly directed flow along its boundary (Steiner et al. 1998). The presence of shock fronts could also imply discontinuities in the velocity field. All these sudden variations of the physical quantities describing the atmosphere might reflect into sudden variations in absorption and emissivity.
Radiative transfer calculations make use of discrete atmospheric models and one can distinguish between two cases: on the one hand, there are semiempirical models (e.g., Vernazza et al. 1981; Fontenla et al. 1990, 1993), which are still frequently used to synthesize spectral lines. In these models, the physical parameters across the atmosphere usually vary smoothly. On the other hand, stateoftheart 3D RMHD simulations of stellar atmospheres show contact discontinuities, steep gradients in state variables, and shock fronts. In fact, if the MHD equations are solved by the method of characteristics then shocks are represented by mathematical discontinuities (Mihalas 1981), while in the case of a nonideal fluid, sharp structures might be smeared over a few zones in the mesh. Nonetheless, one still deals with a quasi discontinuity when the jump in physical states occurs on a spatial scale much smaller than that of the numerical discretization (Steiner et al. 2016).
Discrete models provide the physical parameters describing the atmosphere at a discrete set of spatial points only. Section 2.1 explains that, in this case, there is no well defined concept of discontinuity. However, one is still interested in determining the possible impact of such features on the numerical solution.
4.2. Radiative transfer equation
The monochromatic (timedependent) transfer of unpolarized light is described by the following scalar partial differential equation (Mihalas 1978)
where I_{ν} is the specific intensity of radiation propagating in the ray path direction at time t. Moreover, χ_{ν} and ϵ_{ν} are the absorption and the emission coefficients, respectively, ν is the frequency, and c is the speed of light. The spatial coordinate s denotes the position along the ray under consideration. Equation (11) is a hyperbolic partial differential equation and, specifically, an advection equation with source and sink terms. It is known that computing the numerical solution of advection equations can be particularly challenging. However, in the vast majority of radiative transfer problems in astrophysics, the photon freeflight time is much smaller than the fluid dynamics timescales and the radiation field fully stabilizes before any change occurs in the medium. Therefore, one usually assumes the steadystate version of Eq. (11), consisting in the linear firstorder inhomogeneous ODE given by (Mihalas 1978)
Assuming that F_{ν} is Riemann integrable in the interval [s_{0}, s], one formally integrates Eq. (12) and obtains
Here, the problem of calculating the emergent radiation consists in the evaluation of the integral in the righthand side of Eq. (13), which in turn depends on the specific intensity I_{ν}.
When dealing with radiative transfer in astrophysical plasmas, it is common to refer to the concept of optical depth. Replacing the coordinate s by the optical depth τ_{ν} defined by
one recasts Eq. (12) into
where S_{ν} = ϵ_{ν}/χ_{ν} is the source function, that is, the ratio between the emission and absorption coefficients. Moreover, the formal integration of Eq. (15) in the interval [τ_{ν, 0},τ_{ν}] gives
where τ_{ν, 0} ≥ τ_{ν}. Here, the problem of calculating the emergent radiation reduces to an integral evaluation.
4.3. Numerical approximations
Auer (2003) defined the (numerical) formal solution of Eq. (12) as the evaluation of the radiation field, given knowledge of the opacity χ_{k} and the emissivity ϵ_{k} coefficients at each grid point s_{k}, and of the boundary conditions.
The first step in the numerical approach to an ODE involves the discretization of the integration domain. Therefore, one discretizes the ray path under consideration with a spatial depth grid {s_{k}} (or an optical depth grid {τ_{k}}) with k = 0, …, N. The index k increases along the propagation direction. The numerical approximation of a certain quantity at node s_{k} (or τ_{k}) is indicated by substituting the explicit dependence on s (or τ) with the subscript k, for instance,
For notational simplicity, the frequency dependence is omitted. Inserting numerical quantities in Eq. (13), one obtains^{1}
Different approximations of the integral on the righthand side yield different numerical schemes. Wellknown examples are the backward Euler method, the implicit trapezoidal method, and the thirdorder RungeKutta method (see Table 1). Note that some numerical schemes (e.g., RungeKutta 3) use intermediate grid points for the integration. In this case, one must recover the relevant quantities at offgrid points, which are typically provided through interpolation. This procedure may alter absorption and emission coefficients and introduce numerical errors.
Numerical schemes.
However, the most common formal solutions ground on Eq. (16). Inserting numerical quantities, one gets
where, from Eq. (14), one has
This recursive relation is the starting point for the socalled exponential integrators (Hochbruck & Ostermann 2010), which are a family of numerical schemes that exhibit strong stability properties^{2}. The numerical problem is reduced to the quadrature of the integrals in Eqs. (17) and (18), for which different techniques are available. The knowledge of S and χ at the neighboring grid points suggests the use of interpolatory quadratures (see Sect. 3.3). The standard choice is to use Lagrange polynomials and the simplest strategy is either constant or linear approximations inside the integration interval. However, Auer & Paletou (1994) noted that in order to recover the diffusion approximation, it is necessary to use parabolic or higherorder interpolations of S. Olson & Kunasz (1987) and Kunasz & Olson (1988) made use of parabolic interpolations and, alternatively, Mihalas et al. (1978) proposed cubic Hermite interpolations of S.
Section 2 explains that the accuracy of numerical schemes for ODEs rely on different assumptions on the degree of smoothness of the functional F. Moreover, Sect. 3.2 argued that polynomial interpolations are inaccurate around discontinuities, whereas Sect. 3.3 showed that the accuracy of interpolatory quadratures relies on a sufficient smoothness of the integrand.
4.4. Polarized radiative transfer equation
Most classical radiative transfer problems do not consider polarization. However, the transfer of polarized light is of particular interest in many applications. The radiative transfer of partially polarized light is described by the system of firstorder coupled inhomogeneous ODEs given by
where s is the spatial coordinate measured along the ray under consideration, I_{ν} = (I, Q, U, V)^{T} is the Stokes vector,
is the 4 × 4 propagation matrix, and ϵ_{ν} = (ϵ_{I}, ϵ_{Q}, ϵ_{U}, ϵ_{V})^{T} is the emission vector (Landi Degl’Innocenti & Landolfi 2004). For notational simplicity the frequency dependence of all the vectorial and matricial entries is omitted.
The peculiarity of the polarized radiative transfer equation originates simply from its matricial character (Landi Degl’Innocenti & Landolfi 2004). In fact, the assumption of vanishing offdiagonal coefficients in the propagation matrix decouples Eq. (19) into four independent scalar problems formally identical to Eq. (12), reducing the vectorial problem to the scalar one. The diagonal element of the propagation matrix corresponds to the total absorption coefficient in the unpolarized case, that is, η_{I} ≡ χ_{ν}.
The generalization of the definition of the scalar formal solution to the polarized case consists in substituting radiation intensity, opacity, and emissivity by Stokes vector, propagation matrix, and emission vector, respectively. The polarized problem is considerably more complex than the scalar one. For instance, the problem cannot be reduced to simple quadratures (see Landi Degl’Innocenti & Landolfi 2004). Moreover the radiative transfer equation for polarized light exhibits stiff behavior, and numerical schemes face instability issues (Janett & Paganini 2018). An extensive analysis of the numerical solution of the polarized radiative transfer equation is given by Janett et al. (2017a,b), who characterized several exponential integrators (DELO methods) and highorder formal solvers.
Just like in the scalar case, the vectorial formal solution suffers from numerical issues in the presence of discontinuities. For the sake of simplicity, the full convergence analysis is not repeated, but the numerical tests presented in Sect. 5 comprehend the polarized problem.
4.5. Two and threedimensional problems
In 1D problems, any ray under consideration is discretized by the grid points. However, in 2D or 3D problems, the ray path (long characteristics) or the downwind and upwind ray segments (short characteristics) do not generally coincide with the grid points at which the physical quantities are given. At each integration step, one must then recover the neighboring absorption and emission (or source function) coefficients and the upwind intensity I_{k} at offgrid points^{3}. This is typically performed through interpolations, whose role is to “fill in the gaps” between the known discrete values (Auer 2003; Steffen 2017).
In 2D problems, the usual strategy is to discretize the ray by taking its intersections with the segments connecting the grid points, dealing then with 1D interpolations to recover the offgrid points quantities (Auer 2003). The linear interpolations from the values of the nearest grid points is the simplest and fastest method, but it leads to a large numerical dispersion (Fabiani Bendicho 2003; Ibgui et al. 2013). As the accuracy of the formal solution depends on the accuracy of the interpolated offgrid points quantities, it is often necessary to use highorder interpolations.
In 3D problems, an efficient approach is to discretize the ray by taking its intersections with the individual planes defined by the grid points, dealing then with 2D interpolations. The simplest approach is to use bilinear interpolations, because a cell wall has four corner points. Once again, linear interpolations are not suitable to reach high accuracy. Extending to higherorder, one can use biquadratic or bicubic interpolations (Dullemond 2013).
Section 3.2 explains that if the behavior of the quantities to be interpolated is discontinuous or particularly intermittent, highorder interpolations oscillate and may introduce artifacts. Both Kunasz & Auer (1988) and Auer (2003) highlighted that highorder Lagrange polynomials may introduce spurious extrema that could lead to nonphysical negative values of the source function or of the intensity.
5. Numerical tests
Numerical analysis alone is merely theory. Indeed, numerical evidence is mandatory to assess the mathematical predictions given in the previous sections. Numerical evidence for oscillatory interpolations is easily accessible in the literature and therefore is not exposed here. In order to determine whether and how discontinuities (or highgradients) affect the numerical integration of Eq. (19), different scenarios are investigated. Since Eq. (19) is always integrated along a ray (Auer 2003), the 1D geometry is adopted. At first, different analytical models for the scalar radiative transfer are analyzed and, secondly, the vectorial polarized case is considered. One should avoid drawing premature conclusions on the mere basis of small details of error curves. However, the numerical tests presented in this section allow for some general considerations. Table 1 gives an overview of the numerical schemes used.
5.1. Scalar case
This section presents the integration of Eq. (12) in four different scenarios, where the absorption and emission coefficients show, respectively:

(i)
a smooth behavior;

(ii)
a discontinuity in absorption;

(iii)
a discontinuity in emissivity;

(iv)
a high gradient in absorption.
Spectral line profiles are synthesized using a modified version of the Octave code used by Janett et al. (2017a,b), while the spectral line parameters are identical to those described in Appendix C of Janett et al. (2017a).
Case(i): Smooth atmosphere
In order to verify the predicted order of accuracy of the different numerical schemes, the case of a smooth atmospheric model is first analyzed. Consider Eq. (12) where
Figure 2a shows that the various numerical methods effectively achieve their predicted order of convergence. In this case, the smoothness assumptions are fulfilled and, consequently, the use of highorder methods either allows to reach higher accuracy or makes the use of coarser spatial grids feasible. Moreover, Fig. 2b shows that local errors are bigger where absorption and emission gradients are higher.
Fig. 2. Top row: loglog representation of the global error in the emergent intensity profile as a function of the number of gridpoints. Bottom row: local error as a function of the spatial coordinate s ∈ [0, 100] for a sampling of 50 equispaced grid points for a single wavelength near the line core. The four columns correspond to four different atmospheric models, where the absorption coefficient χ_{ν} and emissivity ϵ_{ν} are given, respectively, by Eqs. (20) and (21) (first column), by Eqs. (22) and (23) (second column), by Eqs. (24) and (25) (third column), and by Eqs. (26) and (27) (fourth column). The parameters are: c_{1} = 0.01, c_{2} = 0.05, c_{3} = 0.15, k_{1} = 1/25, k_{2} = 1/15, k_{3} = 5, j_{1} = 2, j_{2} = 20, s_{d1} = 39.985, and s_{d2} = 39.955. This choice of s_{d1} and s_{d2} warrants a monotonic approach of the grid points to the location of the discontinuity (or sharp gradient) as the grid sampling is refined. The global and the local error are computed as exposed in Eqs. (A.1) (considering the Stokes I component only) and (A.2), respectively, with a reference atmospheric model with 10^{4} grid points. 

Open with DEXTER 
Case(ii): Discontinuity in absorption
The effect of discontinuities in the righthand side of Eq. (12) is first analyzed by considering
A firstorder jump of magnitude j_{1} occurs in the absorption coefficient when the independent variable s reaches the value s_{d1}.
Figure 2c clearly shows the order breakdown due to the firstorder discontinuities. Moreover, Fig. 2d shows that the global error is dominated by the local errors due to the discontinuity.
Case(iii): Discontinuity in emissivity
The effect of discontinuities in the righthand side of Eq. (12) is analyzed by considering
A firstorder jump of magnitude j_{2} occurs in the emission coefficient when the independent variable s reaches the value s_{d1}.
Just like in case (ii), Fig. 2e clearly shows the predicted order breakdown due to the firstorder discontinuities. Moreover, Fig. 2f shows that the global error is dominated by the local errors due to discontinuity.
Case(iv): High gradient in absorption
Atmospheric models can involve rapidly varying physical parameters without an actual mathematical discontinuity. The effect of sharp variations in the righthand side of Eq. (12) is analyzed with the model
A sharp gradient occurs in the absorption coefficient when the independent variable s reaches the value s_{d2}. For k_{3} → ∞ the function tends to a step function with the discontinuity at s = s_{d2}. This case represents the continuous version of case (ii).
For coarse samplings, Fig. 2g clearly shows the order breakdown due to the high gradient. The global error is dominated by the local errors due to the high gradient and such errors present only firstorder convergence until the sampling starts to resolve the sharp gradient. In fact, Sect. 7 explains that there is no difference between a discontinuity and a sharp gradient for a set of discrete data which is not able to resolve the local feature. Figure 2g shows that, once the sampling resolves the high gradients, the numerical methods start to converge according their predicted order of accuracy.
5.2. Polarized case
This section presents some numerical evidence for the integration of Eq. (19) in the case of discontinuities in the atmospheric physical parameters. From these numerical tests, one can still draw conclusions about the scalar problem considering the Stokes I component only.
Stokes profiles of the photospheric Sr I line at 4607.3 Å are synthesized using a modified version of the RH code of Uitenbroek (2001) that allows to switch between different formal solvers and to sequentially perform Stokesprofile syntheses with a set of discrete atmospheric models (see Janett et al. 2018). In these calculations, spectral lines are synthesized under the assumption of LTE conditions and the polarization is produced by the Zeeman effect alone.
Figure 3 shows the considered atmospheric models. These are modified versions of the model C of Fontenla et al. (1993), the socalled FALC model^{4}. These models present, respectively:

(i)
a smooth profile;

(ii)
a discontinuity in the temperature;

(iii)
a discontinuity in the line of sight velocity;

(iv)
a discontinuity in the magnetic field.
Fig. 3. Temperature (top), upward directed line of sight velocity (middle), and three components of the magnetic field (bottom), all as a function of the geometrical height in the atmosphere. The dashed lines represent the FALC smooth atmospheric model, while solid lines show discontinuities in temperature, velocity, and magnetic field profiles. 

Open with DEXTER 
These atmospheric models are homogeneously resampled in logτ_{c} (τ_{c} is the continuum optical depth at the wavelength λ = 5000 Å with different gridpoint densities in the optical depth interval −8.6 ≤ logτ_{c} ≤ 1.4, which encompasses 10 decades. The Sr I line at 4607.3 Å is sampled with around 500 points equispaced in frequency in a spectral interval of a few angstroms around the core.
Figure 4 gives the loglog representation of the global error in the emergent Stokes profiles as a function of the number of pointsperdecade of continuum optical depth, i.e., the number of grid points that sample a variation of one order of magnitude in the optical depth at the wavelength λ = 5000 Å. These error curves allow for different general considerations.
Fig. 4. loglog representation of the global error for the Stokes vector components I, Q, U and V as a function of the number of pointsperdecade of continuum optical depth for different formal solvers (color coded). The Sr I line at 4607.3 Å is considered for the smooth FALC atmospheric model (first row), and for three additional atmospheric models that contain discontinuities in temperature (second row), in line of sight velocity (third row), and in magnetic field (fourth row). The atmospheric models are described in Fig. 3. The global error is computed as exposed in Eq. (A.1). 

Open with DEXTER 
In the preasymptotic regime (below 2 pointsperdecade), the comparison between the rows of Fig. 4 does not reveal essential differences in the error curves. In this regime, the accuracy of numerical schemes strongly depends on the specific sampling.
Asymptotic convergence rates become relevant above 3 pointsperdecade. The first row of Fig. 4 shows that the different numerical methods effectively achieve the predicted (second, third, and fourth) order of accuracy. In this case, the smoothness assumptions are fulfilled, allowing the asymptotic order of accuracy of the numerical solution and highorder schemes outperform loworder methods.
By contrast, the second, third, and fourth rows of Fig. 4 clearly exhibit the order breakdown due to the discontinuity. In fact, all the numerical methods drop to firstorder convergence, making the application of highorder schemes pointless. This attests to the fact that discontinuities in the atmospheric physical parameters effectively induce firstorder discontinuities in Eq. (19) and affect the convergence (and the accuracy) of numerical solutions. In the fourth row, the error curve of Stokes I seems to be only mildly affected by the discontinuity. This is because the magnetic field has a meager impact on Stokes I with respect to Q, U, and V.
Droppedtofirstorder methods require among 30–40 pointsperdecade to achieve an accuracy of E_{i} ≤ 10^{−2}, for i = 1, 2, 3, 4 (see Eq. (A.1)), whereas they require prohibitive numerical grids to produce E_{i} ≤ 10^{−3} accurate spectra. This demonstrates the need for highorder wellbehaved formal solvers able to face discontinuities in the radiative transfer equation.
6. Numerical methods for discontinuous ODEs
The need for efficient and accurate numerical methods for treating discontinuous differential systems has been widely recognized. Most of the literature on the topic is focused on adaptive stepsize numerical methods, that is, numerical schemes that dynamically choose the length of steps along the integration. In this way, the algorithm maintains the desired accuracy of the approximation by decreasing the step size when large derivatives (or singularities) appear and improves efficiency by increasing the step size whenever the smoothness of the system allows it. In reallife radiative transfer problems, the sampling of the atmospheric model is given. This implies that the step size is fixed a priori and thus the adaptive stepsize technique cannot be applied. In principle, one might circumvent this problem by providing the absorption and emission quantities at intermediate gridpoints through interpolations. However, Sect. 3 points out that interpolations suffer from local oscillations around discontinuities.
When facing a discontinuous ODE, Gear & Osterby (1984) distinguish among four main stages in a numerical procedure:

detecting the discontinuity,

locating the discontinuity,

crossing the discontinuity,

restarting after the discontinuity.
Moreover, Dieci & Lopez (2012) classified the different numerical schemes in three main categories: time stepping methods, event driven methods, and regularization of the differential system. This section briefly reviews the numerical schemes for the treatment of discontinuous ODEs, investigating their applicability to Eqs. (12) and (19) in discontinuous media.
6.1. Time stepping methods
Time stepping methods usually rely on local error estimators to detect and locate discontinuities: a large value of the local error estimate indicates the presence of a discontinuity and different strategies to estimate local errors are available (Shampine & Watts 1971; Butcher & Johnston 1993). Alternatively, Calvo et al. (2003) exposed the possibility to detect discontinuities by using the socalled pairs of embedded forms or by monitoring the size of an estimate of the local defect. Some discontinuities might avoid detection. These would probably have negligible effects on the local truncation error and they can consequently be ignored. In contrast, for some smooth but rapidly varying problems, one may detect discontinuities that mathematically do not exist. In this case, it is often reasonable to assume that a discontinuity is present and to act accordingly (Gear & Osterby 1984).
In radiative transfer problems, one can determine at each step whether or not the numerical solution is crossing a discontinuity. This detection should be computationally inexpensive and must fit well into the code, that is, it should be based on quantities that the code computes at each step. Alternatively, one may locate discontinuities using the discrete values of absorption and emission coefficients and by detecting the discontinuities in the discrete data before the integration or during the integration. This approach involves data analysis along the entire ray path and has the advantage that the treatment of the discontinuity is performed in a separate preliminary step. A more indepth study is required to support this strategy, but some suitable techniques are already available (Gelb & Tadmor 1999, 2002; Mallat & Hwang 2006; Oeffner et al. 2013).
In reallife radiative transfer problems, the location of a discontinuity is, in some sense, harder to decipher than commonly assumed. In discrete data, it is usually impossible to exactly locate the discontinuity, because there is a lack of knowledge and the jump can occur anywhere inside the interval.
After having located a discontinuity, time stepping methods adapt time steps to ensure that the local error is kept below a specified tolerance ϵ. This automatically enforces accuracy in spite of the decreased smoothness of the problem. Accordingly, Dieci & Lopez (2012) defined the passing stepsize Δt̃ as
where K_{1} is the firstorder jump in f defined in Sect. 2.2. However, radiative transfer problems usually assume a priori fixed step sizes and the adaptive techniques are not suitable.
After having crossed a discontinuity in [t_{k − 1}, t_{k}], the numerical methods start up the usual integration. However, one must be careful about the previously obtained information, because back values of y and f (i.e., y_{k − 1}, f_{k − 1}, y_{k − 2}, f_{k − 2}, …) may not correspond to those of smooth functions. In case of a discontinuity of order q, they differ by terms of the order 𝒪(h^{q}) and 𝒪(h^{q−1}), respectively. The use of back values in the numerical integration give rise to terms in the local truncation error of order 𝒪(h^{q}), and clearly, it is always safe to use a numerical method that does not use them (e.g., the trapezoidal method and DELOlinear).
6.2. Eventdriven methods
In many problems, discontinuities do not occur completely unexpectedly. The socalled eventdriven methods assume that there is an event function, that is, a function that changes sign at the location of the discontinuity (Dieci & Lopez 2012). Mao & Petzold (2002) pointed out that most discontinuitydetection algorithms are based on the rootfinding of event functions. In fact, when the numerical solution reaches a root of the event function the discontinuity is detected and located. Once the location of the discontinuity is precisely known, then the numerical method includes it as an extra grid point and the algorithm accordingly adjusts the stepping and restarts at that point. In doing so, the local smoothness assumptions are fulfilled, maintaining the asymptotic correctness of the numerical solution and avoiding the order breakdown (Mannshardt 1978).
Unfortunately, no event function is known in the radiative transfer problems and eventdriven methods are consequently not suitable in this context.
6.3. Regularization of the differential system
An alternative strategy is to regularize (or smoothen) the differential system, removing the mathematical discontinuities. Undoubtedly, this leads to simplifications in the theory on ODE^{5}. A more technical and thorough investigation of this topic can be found in Llibre & Teixeira (1997) and Llibre et al. (2006). However, due to the large derivatives that replace the structural discontinuities, the regularized system becomes quite stiff and small time steps are then required. Moreover, regularization can lead to changes in the dynamics of the original nonsmooth system. For these reasons, the regularization of the differential system strategy is not particularly adequate for radiative transfer problems.
7. Interpolations for discontinuous discrete data
Different approximations can be used to recover the various intermittent quantities involved in radiative transfer problems at offgrid points. The goal is to obtain highorder accuracy in smooth regions, avoiding spurious effects around sharp gradients and discontinuities. The standard strategy is to use a locally monotonic^{6} interpolation for each direction (x, y, z and the shortcharacteristic direction in Cartesian grids).
This section shortly reviews different interpolation (or reconstruction) techniques for scattered quantities, providing some illustrative applications to the radiative transfer problem conducted in the past. In particular: cubic Hermite splines, Bézier curves, piecewise rational polynomials, slope limiters, and highorder (weighted) essentially nonoscillatory approximations.
7.1. Cubic Hermite splines
Given a set of points {x_{i}}, the Hermite interpolation H does not only match a set of function values {y_{i}}, but also its derivatives (see Janett et al. 2017b, Sect. 4). For notational simplicity one defines the normalized variable t ∈ [0, 1] as
The cubic Hermite interpolation, approximating the function y(t) inside the interval t ∈ [0, 1], reads
It is known that the cubic Hermite interpolant is fourthorder accurate if the derivatives are at least thirdorder, and thirdorder accurate if the derivatives are secondorder, and so on (Dougherty et al. 1989). Therefore, assuming derivatives that are at least thirdorder accurate, the error scales as 𝒪(h^{4}), indicating that the cubic Hermitian interpolation is fourthorder accurate. Monotonicity can be ensured by a suitable choice of the derivatives. The main weakness of this approach is that it necessarily degenerates to a linear interpolation near smooth extrema (Shu 1998).
Fritsch & Carlson (1980) proposed a twopass algorithm to recover such derivatives. Later on, Fritsch & Butland (1984) and Steffen (1990) added alternative formulas to recover secondorder accurate first derivatives^{7}. Auer (2003) suggested the use of monotonic Hermite interpolants in radiative transfer problems and Ibgui et al. (2013) used monotonic cubic Hermite polynomials (version of Fritsch & Butland 1984) in the IRIS code.
7.2. Bézier curves
These wellknown interpolations make use of the socalled control points (or weights) to suppress spurious extrema. A Bézier curve of degree q applied to a set of function values {y_{i}} at positions {x_{i}} inside the interval t ∈ [0, 1] can be defined as
where C_{n} are the control points, and the Bernstein polynomials B_{n, q}(x) are given by
The first and the last control points define the start and end points of the Bézier curve in the interval, that is,
All the remaining points, conventionally referred to as weights, are usually used to shape the curve. Moreover, a Bézier curve always lies in the convex hull of the control points, that is, in the smallest set that contains the line segment joining every pair of control points (see Janett et al. 2017b, Sect. 5). One can avoid the creation of new extrema by adjusting the weights. However, the highorder accuracy of Bézier interpolations is achieved by forcing the Bézier interpolants to be identical to the corresponding degree Hermite interpolants but, in this case, monotonicity is not guaranteed.
Auer (2003) suggested the use of Bézier curves in radiative transfer problems, because of their suitability for preventing spurious behavior near rapid changes in the absorption and emission coefficients. Štěpán & Trujillo Bueno (2013) used Bézier interpolations in the PORTA code, whereas de la Cruz Rodríguez & Piskunov (2013) made use of them to construct the quadratic and cubic DELOBézier formal solvers.
7.3. Piecewise rational polynomials
As an alternative to the use of polynomials for the interpolation of monotonic data, Gregory & Delbourgo (1982) and Delbourgo & Gregory (1983) opted for the application of piecewise rational quadratic functions. Delbourgo & Gregory (1985) constructed a monotone, piecewise rational cubic interpolation, that includes the rational quadratic function as a special case.
A piecewise rational cubic function R(t), approximating the function y(t) with t ∈ [0, 1], reads
where
and
The condition r_{i} > −1 ensures a strictly positive denominator in the rational cubic, while r_{i} = 3 reduces the rational cubic to the standard cubic Hermite polynomial given in Eq. (28). Here, the parameter r_{i} is chosen to ensure that the interpolant preserves monotonicity. In particular, if
then the rational cubic defined by Eq. (29) reduces to the rational quadratic form
which yields a monotonic interpolant.
In most applications, the derivatives and are not known in advance and hence must be numerically determined. Delbourgo & Gregory (1985) proposed different approximations for the derivative parameters and showed that an error 𝒪(h^{4}) can be expected when 𝒪(h^{3}) derivative information is given at the data points. Finally, the rational quadratic given by Eq. (32) can be used to construct a C^{2} rational spline which interpolates strictly monotonic data (Delbourgo & Gregory 1985).
7.4. Slope limiters
An alternative method to avoid spurious oscillations near discontinuities is to reduce the order of accuracy of the interpolation, for example using a linear rather than a quadratic interpolant near the discontinuity, or by reducing the slope of reconstructions applying limiters. In fact, large gradients and discontinuities produce local extrema in highorder reconstructions, inducing a lack of monotonicity (overshoots), demonstrating the necessity to limit the spatial derivatives inside each cell to physically meaningful values. This is already guaranteed for the Godunov (constant) scheme, but some limitations are required for the Van Leer (linear) and the parabolic reconstructions. In computational fluid dynamics, this problem is solved by the socalled slope limiters. There exists a variety of slope limiters and broad literature attests to the importance of this topic (e.g., Colella & Woodward 1984; Sweby 1984; Berger et al. 2005; Colella & Sekora 2008; Velechovský et al. 2013). One disadvantage of this approach is that it necessarily degenerates to a linear interpolation near smooth extrema (Shu 1998).
Steiner et al. (2016) introduced the use of limiters for the interpolations and reconstructions of the coefficients involved in the polarized radiative transfer equation, purposely taking discontinuities into account.
7.5. Essentially nonoscillatory and weighted essentially nonoscillatory approximations
Provided the smoothness of the function inside the interval, it is known that the wider the stencil, i.e., the set of grid points considered, the higher the order of accuracy of the interpolation. The most common interpolations (see Sects. 7.1–7.3) are based on fixed stencils. However, fixed stencil interpolations of second or higher order are oscillatory in the presence of a discontinuity.
Essentially nonoscillatory (ENO) methods are based on a nonlinear adaptive procedure that automatically chooses the stencil that leads to the locally smoothest interpolant, avoiding crossing discontinuities. A broad range of literature on the topic is available (Liu et al. 1994; Shu 1998, 2009). The classical ENO scheme chooses, among several candidates, the smoothest stencil to work with and discards the rest. This stencil is then used to construct a polynomial interpolation. Otherwise, weighted ENO (WENO) schemes use weighted linear combinations of smallerstencil polynomials to obtain higherorder accuracy. WENO approximations guarantee a nonoscillatory result since the contribution from any stencil containing the discontinuity has an essentially zero weight.
Both ENO and WENO schemes are especially suitable for problems containing both discontinuities and smooth structures. Moreover, they do not necessarily degenerate to a linear interpolation near smooth extrema.
As an explicit example, an illustrative thirdorder accurate ENO interpolation is presented in the following.
To obtain a thirdorderaccurate interpolation of the function y(x), one needs a threepoint stencil. Thus in cell [x_{i}, x_{i + 1}], one starts with the twopoint stencil S_{2} = {x_{i}, x_{i + 1}}. Subsequently, one has two choices to expand the stencil: by adding either the left neighbor x_{i − 1} or the right neighbor x_{i + 2}. The selection criteria is to compare the local smoothness of the function to be interpolated, measured in terms of divided differences (Shu 1998). For instance, one takes the decision by comparing the absolute values of the two relevant divided differences y[x_{i − 1}, x_{i}, x_{i + 1}] and y[x_{i}, x_{i + 1}, x_{i + 2}]. A smaller one implies that the function is smoother in that stencil. Therefore, if
the threepoint stencil is taken as
otherwise,
Once the interpolation stencil is determined, one proceeds with a standard (Lagrange) polynomial interpolation. Clearly, one can continue this procedure to add more grid points to the stencil, constructing higherorder ENO approximations.
A more thorough investigation on the possible applications of ENO and WENO techniques in radiative transfer problems is in progress.
8. Conclusions
This paper identifies and exposes the relevant problems appearing in the numerical treatment of the radiative transfer equation in discontinuous media. The main aim is not to provide the ultimate numerical method able to face discontinuous problems, but to better understand the specific situations where discontinuity issues appear.
The first part pays particular attention to assumptions and limitations of the standard convergence analysis of numerical schemes for ODEs and interpolations. In fact, standard numerical methods (and interpolations) rely on smoothness assumptions on the functions to be integrated (interpolated) and may perform very inefficiently in the presence of discontinuities, which may drastically increase local errors, reducing the accuracy of the solution and thwarting highorder convergence.
The numerical tests (performed for both the scalar and the polarized case) corroborate analytical predictions: discontinuities in absorption and emission coefficients effectively affect the convergence and, consequently, the accuracy of numerical solutions, inducing the order breakdown where all the numerical methods drop to firstorder convergence. Moreover, local errors are higher around highgradients and unresolved highgradients behave like discontinuities. It is also shown that discontinuities in the atmospheric physical parameters effectively induce firstorder discontinuities in the radiative transfer equation, making the application of highorder schemes pointless. Unfortunately, firstorder convergence is not sufficiently accurate in many situations. In the presence of a discontinuity, the order breakdown in the formal solution is inevitable and the only practical solution is to increase resolution in the atmospheric model. Discontinuities do not cause numerical instability (in the sense of magnification of errors) in the formal solution, but they just decrease accuracy.
The final part summarizes the existing numerical methods for the treatment of discontinuous ODEs and the interpolation techniques for discontinuous discrete data. Standard numerical schemes for discontinuous ODE are shown to be unsuitable for common RTE problems. By contrast, many suitable interpolation techniques are already available. In particular, Essentially NonOscillatory (ENO) and Weighted ENO (WENO) techniques are highorder robust approximations that resolve discontinuities in an accurate and nonoscillatory fashion. An investigation into their possible applications to radiative transfer problems is in progress.
The very same procedure can apply starting from Eq. (15) instead of Eq. (12). Janett & Paganini (2018) analyze the differences between the two cases.
Positiveness of Δτ_{k} guarantees Lstability (Janett & Paganini 2018).
The formula by Fritsch & Butland (1984) is secondorder accurate on uniform grids only and it drops to firstorder on nonuniform grids.
Acknowledgments
The financial support by the Swiss National Science Foundation (SNSF) through grant ID 200021_159206 is gratefully acknowledged. Special thanks are extended to F. Calvo and C. Bertoni for particularly enriching discussions, and to E. Alsina Ballester, L. Belluzzi, and O. Steiner for reading and commenting on previous versions of the paper.
References
 Acary, V., & Brogliato, B. 2008, Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics (Springer Verlag), Lecture Notes Appl. Comput. Mech., 35, 526 [Google Scholar]
 Auer, L. 2003, in Stellar Atmosphere Modeling, eds. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 3 [NASA ADS] [Google Scholar]
 Auer, L. H., & Paletou, F. 1994, A&A, 285, 675 [NASA ADS] [Google Scholar]
 Berger, M., Aftosmis, M. J., & Murman, S. M. 2005, Analysis of Slope Limiters on Irregular Grids 43rd AIAA Aerospace Sciences Meeting, 2361 [Google Scholar]
 Butcher, J., & Johnston, P. 1993, J. Comput. Appl. Math., 45, 203 [Google Scholar]
 Calvo, M., Montijano, J. I., & Rández, L. 2003, Numer. Algorithms, 33, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Calvo, M., Montijano, J., & Rández, L. 2008, Bol. Soc. Esp. Math. Apl., 44, 33 [Google Scholar]
 Chartres, B., & Stepleman, R. 1972, SIAM J. Numer. Anal., 9, 476 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Sekora, M. D. 2008, J. Comput. Phys., 227, 7069 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 de la Cruz Rodríguez, J., & Piskunov, N. 2013, ApJ, 764, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Delbourgo, R., & Gregory, J. A. 1983, IMA J. Numer. Anal., 3 [Google Scholar]
 Delbourgo, R., & Gregory, J. A. 1985, SIAM J. Sci. Stat. Comput., 6, 967 [Google Scholar]
 del Toro Iniesta, J. C., & Ruiz Cobo, B. 2016, Sol. Phys., 13, 4 [Google Scholar]
 Dieci, L., & Lopez, L. 2012, J. Comput. Appl. Math., 236, 3967 [Google Scholar]
 Dougherty, R., Edelman, A., & Hyman, J. M. 1989, Math. Comput., 52, 471 [Google Scholar]
 Dullemond, C. P. 2013, Radiative Transfer in Astrophysics. Theory, Numerical Methods and Applications, Tech. Rep. (University of Heidelberg) [Google Scholar]
 Fabiani Bendicho, P. 2003, in Stellar Atmosphere Modeling, eds. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 419 [NASA ADS] [Google Scholar]
 Fontenla, J. M., Avrett, E. H., & Loeser, R. 1990, ApJ, 355, 700 [NASA ADS] [CrossRef] [Google Scholar]
 Fontenla, J. M., Avrett, E. H., & Loeser, R. 1993, ApJ, 406, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Fritsch, F. N., & Butland, J. 1984, SIAM J. Sci. Stat. Comput., 5, 300 [Google Scholar]
 Fritsch, F. N., & Carlson, R. E. 1980, SIAM J. Numer. Anal., 17, 238 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Gear, C. W., & Osterby, O. 1984, ACM Trans. Math. Softw., 10, 23 [Google Scholar]
 Gelb, A., & Tadmor, E. 1999, Appl. Comput. Harmon. Anal., 7, 101 [Google Scholar]
 Gelb, A., & Tadmor, E. 2002, ESAIM Math. Model. Numer. Anal., 36, 155 [Google Scholar]
 Gouzé, J.L., & Sari, T. 2002, Dyn. Syst., 17, 299 [Google Scholar]
 Gregory, J. A., & Delbourgo, R. 1982, IMA J. Numer. Anal., 2, 123 [Google Scholar]
 Hochbruck, M., & Ostermann, A. 2010, Acta Numer., 19, 209 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013, A&A, 549, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Janett, G., & Paganini, A. 2018, ApJ, 857, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Janett, G., Carlin, E. S., Steiner, O., & Belluzzi, L. 2017a, ApJ, 840, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Janett, G., Steiner, O., & Belluzzi, L. 2017b, ApJ, 845, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Janett, G., Steiner, O., & Belluzzi, L. 2018, ApJ, 865, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Kunasz, P., & Auer, L. H. 1988, J. Quant. Spectr. Rad. Transf., 39, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Kunasz, P. B., & Olson, G. L. 1988, J. Quant. Spectr. Rad. Transf., 39, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Landi Degl’Innocenti, E., & Landolfi, M. 2004, in Polarization in Spectral Lines (Dordrecht: Kluwer Academic Publishers), Astrophys. Space Sci. Lib., 307 [Google Scholar]
 Landry, C., Caboussat, A., & Hairer, E. 2009, SIAM J. Sci. Comput., 31, 3806 [Google Scholar]
 Liu, X.D., Osher, S., & Chan, T. 1994, J. Comput. Phys., 115, 200 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Llibre, J., & Teixeira, M. A. 1997, Dyn. Syst., 3, 235 [Google Scholar]
 Llibre, J., Silva, P. R., & Teixeira, M. A. 2006, J. Dyn. Diff. Equ., 19, 309 [Google Scholar]
 Louis, R. E., Bellot Rubio, L. R., Mathew, S. K., & Venkatakrishnan, P. 2009, ApJ, 704, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Mallat, S., & Hwang, W. L. 2006, IEEE Trans. Inf. Theor., 38, 617 [Google Scholar]
 Malmborg, J., & Bernhardsson, B. 1997, Nonlinear Anal., 30, 337 [Google Scholar]
 Mannshardt, R. 1978, Numer. Math., 31, 131 [Google Scholar]
 Mao, G., & Petzold, L. 2002, Comput. Math. Appl., 43, 65 [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (San Francisco: W.H. Freeman and Company) [Google Scholar]
 Mihalas, D. 1981, J. Quant. Spectr. Rad. Transf., 25, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., Auer, L. H., & Mihalas, B. R. 1978, ApJ, 220, 1001 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, D. A. N., Schlichenmaier, R., Steiner, O., & Stix, M. 2002, A&A, 393, 305 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oeffner, P., Sonar, T., & Wirz, M. 2013, Appl. Math., 4, 1 [Google Scholar]
 Olson, G. L., & Kunasz, P. B. 1987, J. Quant. Spectr. Rad. Transf., 38, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Richards, F. 1991, J. Approx. Theory, 66, 334 [Google Scholar]
 Shampine, L. F., & Watts, H. A. 1971, Math. Comput., 25, 445 [Google Scholar]
 Shu, C.W. 1998, in Essentially Nonoscillatory and Weighted Essentially Nonoscillatory Schemes for Hyperbolic Conservation Laws, ed. A. Quarteroni (Berlin, Heidelberg: Springer Berlin Heidelberg), 325 [Google Scholar]
 Shu, C.W. 2009, SIAM Rev., 51, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Steffen, M. 1990, A&A, 239, 443 [NASA ADS] [Google Scholar]
 Steffen, M. 2017, Mem. Soc. Astron. It., 88, 22 [NASA ADS] [Google Scholar]
 Steiner, O. 2000, Sol. Phys., 196, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Steiner, O., GrossmannDoerth, U., Knölker, M., & Schüssler, M. 1998, ApJ, 495, 468 [NASA ADS] [CrossRef] [Google Scholar]
 Steiner, O., Züger, F., & Belluzzi, L. 2016, A&A, 586, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Štěpán, J., & Trujillo Bueno, J. 2013, A&A, 557, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sweby, P. K. 1984, SIAM J. Numer. Anal., 21, 995 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Tscharnuter, W. 1977, A&A, 57, 279 [NASA ADS] [Google Scholar]
 Uitenbroek, H. 2001, ApJ, 557, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Velechovský, J., Liska, R., & Shashkov, M. 2013, Comput. Fluids, 83, 164 [Google Scholar]
 Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, Z., & Martin, C. F. 1997, J. Comput. Appl. Math., 87, 359 [Google Scholar]
Appendix A: Error calculation
Denoting with I^{ref}(ν) and I^{num}(ν) the reference and the numerically computed emergent Stokes vectors, respectively, at the frequency ν, the global error for each Stokes vector component is computed as
where i = 1, 2, 3, 4 indicates the Stokes parameters (I, Q, U, V). The error is given by the maximal discrepancy between the reference and the simulated Stokes parameter over the spectral interval considered, normalized by the maximal amplitude in the reference profile. The reference emergent Stokes profile is calculated with the cubic Hermitian method (and crosschecked with the cubic DELOBézier method) using a hyperfine grid sampling with 10^{3} pointsperdecade of continuum optical depth.
Denoting with and the reference and the numerically computed specific intensities at the frequency ν and at the grid location s_{k}, the local error is computed as
All Tables
All Figures
Fig. 1. The hyperbolic tangent (top panel, continuous blue) and the modified Heaviside (bottom panel, continuous blue) functions are sampled in the interval x ∈ [−1, 1] with 10 (red circles) and 40 (green dots) grid points. No difference between the coarser (red) samplings is visible, while this is not true for the finer (green) samplings. 

Open with DEXTER  
In the text 
Fig. 2. Top row: loglog representation of the global error in the emergent intensity profile as a function of the number of gridpoints. Bottom row: local error as a function of the spatial coordinate s ∈ [0, 100] for a sampling of 50 equispaced grid points for a single wavelength near the line core. The four columns correspond to four different atmospheric models, where the absorption coefficient χ_{ν} and emissivity ϵ_{ν} are given, respectively, by Eqs. (20) and (21) (first column), by Eqs. (22) and (23) (second column), by Eqs. (24) and (25) (third column), and by Eqs. (26) and (27) (fourth column). The parameters are: c_{1} = 0.01, c_{2} = 0.05, c_{3} = 0.15, k_{1} = 1/25, k_{2} = 1/15, k_{3} = 5, j_{1} = 2, j_{2} = 20, s_{d1} = 39.985, and s_{d2} = 39.955. This choice of s_{d1} and s_{d2} warrants a monotonic approach of the grid points to the location of the discontinuity (or sharp gradient) as the grid sampling is refined. The global and the local error are computed as exposed in Eqs. (A.1) (considering the Stokes I component only) and (A.2), respectively, with a reference atmospheric model with 10^{4} grid points. 

Open with DEXTER  
In the text 
Fig. 3. Temperature (top), upward directed line of sight velocity (middle), and three components of the magnetic field (bottom), all as a function of the geometrical height in the atmosphere. The dashed lines represent the FALC smooth atmospheric model, while solid lines show discontinuities in temperature, velocity, and magnetic field profiles. 

Open with DEXTER  
In the text 
Fig. 4. loglog representation of the global error for the Stokes vector components I, Q, U and V as a function of the number of pointsperdecade of continuum optical depth for different formal solvers (color coded). The Sr I line at 4607.3 Å is considered for the smooth FALC atmospheric model (first row), and for three additional atmospheric models that contain discontinuities in temperature (second row), in line of sight velocity (third row), and in magnetic field (fourth row). The atmospheric models are described in Fig. 3. The global error is computed as exposed in Eq. (A.1). 

Open with DEXTER  
In the text 