Formal Solutions for Polarized Radiative Transfer. I. The DELO Family

A spatial grid $\{{s}_{k}\}$ $(k=0,\ldots ,N)$ is introduced, discretizing the ray path. The spatial coordinate s and the index k increase along the propagation direction. For a given grid point s_k, the points ${s}_{k-1}$ and ${s}_{k+1}$ represent the upwind and downwind points, respectively, guaranteeing ${s}_{k-1}\leqslant {s}_{k}\leqslant {s}_{k+1}$. Applying the fundamental theorem of calculus to Equation (1) in the interval $[{s}_{k},{s}_{k+1}]$, one obtains

$$\begin{eqnarray}&&{\boldsymbol{I}}({s}_{k+1})-{\boldsymbol{I}}({s}_{k})={\int }_{{s}_{k}}^{{s}_{k+1}}{\boldsymbol{F}}(s,{\boldsymbol{I}}(s)){\rm{d}}s.\end{eqnarray} \tag{ 2 }$$

The integral can be approximated by different numerical quadratures and one possibility is to use the trapezoidal rule, i.e.,

$$\begin{eqnarray*}&&{\int }_{a}^{b}f(x){\rm{d}}x\approx \displaystyle \frac{b-a}{2}[f(a)+f(b)].\end{eqnarray*}$$

Approximating the integral in Equation (2) in terms of the trapezoidal rule, one recovers

$$\begin{eqnarray}&&{\boldsymbol{I}}({s}_{k+1})-{\boldsymbol{I}}({s}_{k})\approx \displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}[{\boldsymbol{F}}({s}_{k},{\boldsymbol{I}}({s}_{k}))+{\boldsymbol{F}}({s}_{k+1},{\boldsymbol{I}}({s}_{k+1}))],\end{eqnarray} \tag{ 3 }$$

where ${\rm{\Delta }}{s}_{k}={s}_{k+1}-{s}_{k}$ is the cell width. The numerical approximation of a certain quantity at node s_k is indicated by substitution of the explicit dependence on s with the subscript k, for instance,

$$\begin{eqnarray*}&&{{\boldsymbol{I}}}_{k}\approx {\boldsymbol{I}}({s}_{k}).\end{eqnarray*}$$

Inserting numerical quantities for ${\boldsymbol{I}}$, ${\bf{K}}$, and ${\boldsymbol{\epsilon }}$ in Equation (3) and applying some algebra, one recovers the implicit linear system

$$\begin{eqnarray}&&{{\boldsymbol{\Phi }}}_{k+1}{{\boldsymbol{I}}}_{k+1}={{\boldsymbol{\Phi }}}_{k}{{\boldsymbol{I}}}_{k}+{{\boldsymbol{\Psi }}}_{k+1}+{{\boldsymbol{\Psi }}}_{k},\end{eqnarray} \tag{ 4 }$$

where the matrices ${{\boldsymbol{\Phi }}}_{k}$ and ${{\boldsymbol{\Phi }}}_{k+1}$ and the vectors ${{\boldsymbol{\Psi }}}_{k}$ and ${{\boldsymbol{\Psi }}}_{k+1}$ are given by

$$\begin{eqnarray*}{{\boldsymbol{\Phi }}}_{k} & = & {\bf{1}}-\displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}{{\bf{K}}}_{k},\\ {{\boldsymbol{\Phi }}}_{k+1} & = & {\bf{1}}+\displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}{{\bf{K}}}_{k+1},\\ {{\boldsymbol{\Psi }}}_{k} & = & \displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}{{\boldsymbol{\epsilon }}}_{k},\\ {{\boldsymbol{\Psi }}}_{k+1} & = & \displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}{{\boldsymbol{\epsilon }}}_{k+1}.\end{eqnarray*}$$

Given the upwind Stokes vector ${{\boldsymbol{I}}}_{k}$ at node s_k, one solves the implicit linear system given by Equation (4), finding the emergent Stokes vector ${{\boldsymbol{I}}}_{k+1}$ at node ${s}_{k+1}$.

The trapezoidal method is therefore classified as an implicit method and belongs to both the famous classes of the Runge–Kutta methods and the linear multistep methods.

In order to recover a good approximation, one tries to maintain as small an error as possible. When discussing numerical schemes, one usually refers to two different kinds of errors: the local truncation error and the global error.

The local truncation error is the error introduced by the numerical scheme in a single step, assuming the exact solution at the precedent step. Considering a scalar initial value problem (IVP) of the general form

$$\begin{eqnarray*}y^{\prime} (t) & = & f(t,y(t)),\\ y({t}_{0}) & = & {y}_{0},\end{eqnarray*}$$

supplied by the discrete grid $\{{t}_{k}\}$ $(k=0,\ldots ,N)$, the local truncation error is defined by

$$\begin{eqnarray*}&&{L}_{k}=\parallel {y}_{k}-y({t}_{k})\parallel ,\ \mathrm{with}\ {y}_{k-1}=y({t}_{k-1}),\end{eqnarray*}$$

where $y({t}_{k-1})$ and $y({t}_{k})$ are the exact solutions at nodes ${t}_{k-1}$ and t_k, respectively, and y_k represents the numerical approximation at node t_k after performing the single step from ${t}_{k-1}$ to t_k. The operator $\parallel \cdot \parallel$ represents any suitable vector norm. A numerical method for an ordinary differential equation has order of accuracy p if the local truncation error, which usually depends on the step size denoted by ${\rm{\Delta }}t$, satisfies

$$\begin{eqnarray}&&{L}_{k}\approx O({\rm{\Delta }}{t}^{p+1}),\,\mathrm{with}\ p\geqslant 1.\end{eqnarray} \tag{ 5 }$$

Hence, the larger the order of accuracy, the faster the error is reduced as ${\rm{\Delta }}t$ decreases.

By contrast, the global error is defined as

$$\begin{eqnarray*}&&{E}_{N}=\parallel {y}_{N}-y({t}_{N})\parallel ,\end{eqnarray*}$$

and represents the accumulation of the local truncation error over all the N steps. If the same amount of error is produced at each step, i.e., without any amplification of errors over subsequent steps, then a local truncation error of $O({\rm{\Delta }}{t}^{p+1})$ implies a global error of $O({\rm{\Delta }}{t}^{p})$. Explicitly

$$\begin{eqnarray}&&{E}_{N}\approx \displaystyle \sum _{k=1}^{N}{L}_{k}\approx C\cdot {\rm{\Delta }}{t}^{p}+O({\rm{\Delta }}{t}^{p+1}),\ \mathrm{with}\ p\geqslant 1,\end{eqnarray} \tag{ 6 }$$

using Equation (5) and $N\propto 1/{\rm{\Delta }}t$. The constant C typically depends on the exact solution and on other parameters of the numerical scheme, but is strictly independent of ${\rm{\Delta }}t$. The concept of amplification of errors is related to the idea of stability, which will be shortly discussed. The intuitive connection between local truncation error and global error given by Equation (6) can be generalized to any suitable discrete grid (e.g., Deuflhard & Bornemann 2002), including logarithmically spaced grids. In the case of non-uniform discrete grids, the step size ${\rm{\Delta }}t$ must be replaced by the maximal step size given by

$$\begin{eqnarray*}&&{\rm{\Delta }}\hat{t}=\mathop{\max }\limits_{k=0,\ldots ,N-1}{\rm{\Delta }}{t}_{k},\end{eqnarray*}$$

with ${\rm{\Delta }}{t}_{k}={t}_{k+1}-{t}_{k}$. The expression "suitable discrete grid" indicates that the maximal step size is inversely proportional to the total number of grid points, i.e., ${\rm{\Delta }}\hat{t}\propto {N}^{-1}$.

The power law presented in Equation (6) implies a linear relationship between the logarithms of the global error E_N and the step size ${\rm{\Delta }}t$, i.e.,

$$\begin{eqnarray*}&&\mathrm{log}({E}_{N})\approx p\cdot \mathrm{log}({\rm{\Delta }}t)+\tilde{C},\end{eqnarray*}$$

where $\tilde{C}=\mathrm{log}(C)$. This relation can be appreciated in a log–log plot, where the resulting straight line, often called the signature or error curve, should show a slope equal to p. The trapezoidal method is well known to be second-order accurate (e.g., Frank 2012) and an explicit example, showing its signature, will be presented later in this paper.

Note, however, that definitions (5) and (6) assume sufficiently small ${\rm{\Delta }}t$ in order to avoid pre-asymptotic behavior, i.e., the fact that for large step sizes the data points are more scattered and do not necessarily follow the power law. For instance, a global error of the form ${E}_{N}\approx {C}_{1}\cdot {\rm{\Delta }}t+{C}_{2}\cdot {\rm{\Delta }}{t}^{2}+{C}_{3}\cdot {\rm{\Delta }}{t}^{3}$ is dominated by the first term for small enough ${\rm{\Delta }}t$, but for large step sizes, higher-order terms may tangibly contribute, depending on their constants C₂ and C₃.

For a numerical scheme, stability means that any numerical error introduced at some stage does not blow up in the subsequent steps of the method. The concept of stability is often related to the concept of stiffness. A differential equation is said to be stiff when some numerical methods have to take an extremely small step size to achieve convergence. Therefore, step-size control is also based on stability requirements, because instabilities lead to a deterioration of accuracy. More details about the concept of stability in the numerical treatment of ordinary differential equations can be found, for instance, in Hackbusch (2014).

There are different ways to determine the stability of a numerical scheme and the stability analysis often depends on the considered class of methods. In the following, the common class of the Runge–Kutta methods is considered. The stability of a numerical method is often deduced through the simple autonomous scalar IVP given by

$$\begin{eqnarray}y^{\prime} (t) & = & \lambda y(t),\\ y(0) & = & {y}_{0},\end{eqnarray} \tag{ 7 }$$

with $\lambda \in {\mathbb{C}}$. The solution $y(t)={y}_{0}{e}^{\lambda t}$ converges to zero as $t\to \infty$ for $\mathrm{Re}(\lambda )\lt 0$. A Runge–Kutta method applied to the IVP (7) can be recast into the form

$$\begin{eqnarray*}&&{y}_{k+1}=\phi (\lambda {\rm{\Delta }}t){y}_{k},\end{eqnarray*}$$

where ϕ is called the stability function. The numerical method is said to be stable if, applied to the IVP (7), it converges to zero for $k\to \infty$ and this condition is equivalent to

$$\begin{eqnarray}&&\parallel \phi (\lambda {\rm{\Delta }}t)\parallel \lt 1.\end{eqnarray} \tag{ 8 }$$

Intuitively, this guarantees that any perturbation in the solution is attenuated with the recursive numerical integration. The stability of a method is therefore related to both the step size ${\rm{\Delta }}t$ and the eigenvalue λ, more precisely to the term $\lambda {\rm{\Delta }}t$. The stability region of a Runge–Kutta method is defined as the set of complex values $\lambda {\rm{\Delta }}t$ for which Equation (8) is satisfied.

The stability analysis for the scalar problem given by Equation (7) can be easily generalized to the linear system of ordinary differential equations given by

$$\begin{eqnarray}{\boldsymbol{y}}^{\prime} (t) & = & {\bf{A}}{\boldsymbol{y}}(t),\\ {\boldsymbol{y}}(0) & = & {{\boldsymbol{y}}}_{0},\end{eqnarray} \tag{ 9 }$$

where the d × d matrix ${\bf{A}}$ has a basis of eigenvectors corresponding to the eigenvalues ${\lambda }^{(1)},\ldots ,{\lambda }^{(d)}$. Equation (9) formally corresponds to the homogeneous version of Equation (1). The emission term is deliberately omitted in the stability analysis because of its locality. In fact, the emission term does not affect the propagation of the information, preventing any contribution to the amplification of errors. As shown in Frank (2012), a Runge–Kutta scheme applied to Equation (9) is stable if and only if it guarantees stability once applied to Equation (7), with λ representing any eigenvalue of ${\bf{A}}$. It is important to mention the fact that the eigenvalues of the propagation operator $-{\bf{K}}$ in Equation (1) have always negative real parts. Therefore, A-stability (see below) is a sufficient condition to avoid any instability problem in the formal solution.

Applying the trapezoidal method to the IVP (7), one recovers the following stability function

$$\begin{eqnarray*}&&{\phi }_{{\rm{T}}}(\lambda {\rm{\Delta }}t)=\displaystyle \frac{1+\lambda {\rm{\Delta }}t/2}{1-\lambda {\rm{\Delta }}t/2}.\end{eqnarray*}$$

The stability region for the trapezoidal method is then given by the condition (8) and it is usually displayed as presented in Figure 1(a). If the stability region of a numerical method contains the whole left-hand side of the complex plane, as in the case of Figure 1(a), then the numerical scheme is said to be A-stable. A broad and exhaustive literature is dedicated to the determination of the specific stability regions for the different numerical methods (Dahlquist 1963; Collatz 1966), but a deep digression would stray from the main aim of this paper.

Zoom In Zoom Out Reset image size

A strong limitation of this simplified stability analysis is the assumption of a constant eigenvalue λ in Equation (7). A less restrictive stability analysis shows that variations of λ along the integration path could affect the stability region of the numerical method. For example, the stability function of the trapezoidal method in the interval $[{t}_{k},{t}_{k+1}]$ would read

$$\begin{eqnarray}&&{\phi }_{{\rm{T}}}({\lambda }_{k},{\lambda }_{k+1},{\rm{\Delta }}{t}_{k})=\displaystyle \frac{1+{\lambda }_{k}{\rm{\Delta }}{t}_{k}/2}{1-{\lambda }_{k+1}{\rm{\Delta }}{t}_{k}/2},\end{eqnarray} \tag{ 10 }$$

where ${\lambda }_{k}$ and ${\lambda }_{k+1}$ are the eigenvalues at the positions t_k and ${t}_{k+1}$, respectively. Thus, the stability region depends on both eigenvalues ${\lambda }_{k}$ and ${\lambda }_{k+1}$, as illustrated in Figure 1. However, the variation of the eigenvalue λ along the integration path strongly depends on the spatial scale. In particular, a conversion from geometrical height to optical depth (see Appendix A) mitigates fluctuations of the eigenvalues of the propagation operator $-{\bf{K}}$, supporting the assumption of a constant eigenvalue λ in the stability analysis.

When designing or choosing a numerical scheme, an important point is the amount of computational time and data storage necessary to execute it (see for instance Goldreich 2008). A suitable parallelization strategy of the radiative transfer problem includes the use of multiple central processor unit (CPU) cores and parallelization via domain decomposition and/or in the frequency domain (Štěpán & Trujillo Bueno 2013). The computational cost may act as an important factor in the choice of the appropriate formal solver, especially for large scale applications in which the repetitive integration of the radiative transfer equation plays a leading role.

The computational cost analysis of numerical schemes must be based on the premise that it depends on the specific coding, programming language, compiler, and computer architecture. Therefore, the interpreted Octave language used in this paper is not suitable to reliably determine time costs. Nevertheless, one can make some objective considerations. First of all, common sense suggests keeping the algorithm as sleek as possible, avoiding any unnecessary superstructure. Second, basic floating-point operations are carried out directly on the CPU, whereas elementary functions are usually emulated on a higher level. Correspondingly, the evaluation of, e.g., an exponential is 10–40 times more expensive than a floating-point multiplication (Schörghofer 2015). A third remark is made on the difference between explicit and implicit schemes. An explicit one-step method calculates the updated numerical value ${y}_{k+1}$ directly from the precedent value y_k, i.e.,

$$\begin{eqnarray*}&&{y}_{k+1}=f({y}_{k}),\end{eqnarray*}$$

while implicit one-step methods find ${y}_{k+1}$ by solving an equation of the type

$$\begin{eqnarray*}&&g({y}_{k},{y}_{k+1})=0,\end{eqnarray*}$$

which results in the additional solution of a 4 × 4 implicit linear system, when considering Equation (1).

Concerning the data storage cost, a simple consideration can be made. In the short characteristic strategy the formal solver integrates step by step the radiative transfer equation along the ray path, using only local atmospheric quantities (Auer 2003). Therefore, the small amount of information retained avoids any data storage problem.

In a very general way, the polynomial ${{\boldsymbol{P}}}_{q}$ in Equation (17) is obtained by the Lagrangian interpolation

$$\begin{eqnarray}&&{\bf{S}}(\tau ,{\boldsymbol{I}}(\tau ))\approx \displaystyle \sum _{i=k-q+1}^{k+1}{{\bf{S}}}_{i}{{\ell }}_{i}(\tau ),\end{eqnarray} \tag{ 20 }$$

where ${{\bf{S}}}_{i}=-{{\bf{K}}}_{i}{{\boldsymbol{I}}}_{i}+{\tilde{{\boldsymbol{\epsilon }}}}_{i}$ and the Lagrange basis polynomials ℓ_i are given by

$$\begin{eqnarray*}&&{{\ell }}_{i}(\tau )=\displaystyle \prod _{\displaystyle \genfrac{}{}{0em}{}{k-q+1\leqslant m\leqslant k+1}{m\ne i}}\displaystyle \frac{\tau -{\tau }_{m}}{{\tau }_{i}-{\tau }_{m}}.\end{eqnarray*}$$

Here k indicates an arbitrary node on the discretized ray path. The integral in Equation (18) can then be solved by parts, yielding an implicit linear system of the form

$$\begin{eqnarray}&&{{\boldsymbol{\Phi }}}_{k+1}{{\boldsymbol{I}}}_{k+1}=\displaystyle \sum _{i=k-q+1}^{k}{{\boldsymbol{\Phi }}}_{i}{{\boldsymbol{I}}}_{i}+\displaystyle \sum _{i=k-q+1}^{k+1}{{\boldsymbol{\Psi }}}_{i},\end{eqnarray} \tag{ 21 }$$

where the different coefficients ${{\boldsymbol{\Phi }}}_{i}$ and ${{\boldsymbol{\Psi }}}_{i}$ depend on the chosen polynomial ${{\boldsymbol{P}}}_{q}$ and on the numerical values ${{\bf{K}}}_{{i}}$ and ${\tilde{{\boldsymbol{\epsilon }}}}_{i}$. Provided the previous Stokes ${{\boldsymbol{I}}}_{i}$ for $i=k-q+1,\ldots ,k$, the linear system (21) can be solved to obtain the numerical approximation ${{\boldsymbol{I}}}_{k+1}$ at ${\tau }_{k+1}$. Therefore, once given the boundary condition I₀ at ${\tau }_{0}$, the recursive application of Equation (21) provides the emergent Stokes vector ${{\boldsymbol{I}}}_{N}$ at the end of the ray path.

As first choice, the effective source function is assumed constant inside the interval $[{\tau }_{k},{\tau }_{k+1}]$, i.e.,

$$\begin{eqnarray*}&&{\bf{S}}(\tau ,{\boldsymbol{I}}(\tau ))\approx {{\bf{S}}}_{k+\tfrac{1}{2}},\end{eqnarray*}$$

and can be approximated by the midpoint rule

$$\begin{eqnarray*}&&{{\bf{S}}}_{k+\tfrac{1}{2}}\approx \displaystyle \frac{{{\bf{S}}}_{k}+{{\bf{S}}}_{k+1}}{2}.\end{eqnarray*}$$

Replacing the polynomial ${{\boldsymbol{P}}}_{q}$ in Equation (18) by the constant approximation described above, one can calculate the integral and after some algebra obtain an implicit linear system formally identical to Equation (4), i.e.,

$$\begin{eqnarray}&&{{\boldsymbol{\Phi }}}_{k+1}{{\boldsymbol{I}}}_{k+1}={{\boldsymbol{\Phi }}}_{k}{{\boldsymbol{I}}}_{k}+{{\boldsymbol{\Psi }}}_{k+1}+{{\boldsymbol{\Psi }}}_{k}.\end{eqnarray} \tag{ 22 }$$

The method described by Equation (22), which might be called DELO-constant, is second-order accurate, as shown in Figure 2. The explicit values of the coefficients ${{\boldsymbol{\Phi }}}_{k}$, ${{\boldsymbol{\Phi }}}_{k+1}$, ${{\boldsymbol{\Psi }}}_{k}$, and ${{\boldsymbol{\Psi }}}_{k+1}$ are provided in Appendix B.

Zoom In Zoom Out Reset image size

The next step is to obtain a relation between ${{\boldsymbol{I}}}_{k}$ and ${{\boldsymbol{I}}}_{k+1}$ when approximating the effective source function by a linear interpolation, i.e., Equation (20) with q = 1,

$$\begin{eqnarray*}&&{\bf{S}}(\tau ,{\boldsymbol{I}}(\tau ))\approx \displaystyle \frac{(\tau -{\tau }_{k+1}){{\bf{S}}}_{k}-(\tau -{\tau }_{k}){{\bf{S}}}_{k+1}}{{\rm{\Delta }}{\tau }_{k}},\end{eqnarray*}$$

for τ in the interval $[{\tau }_{k},{\tau }_{k+1}]$. One proceeds with an analytical integration and after some algebra obtains

$$\begin{eqnarray}&&{{\boldsymbol{\Phi }}}_{k+1}{{\boldsymbol{I}}}_{k+1}={{\boldsymbol{\Phi }}}_{k}{{\boldsymbol{I}}}_{k}+{{\boldsymbol{\Psi }}}_{k+1}+{{\boldsymbol{\Psi }}}_{k},\end{eqnarray} \tag{ 23 }$$

which is an implicit linear system formally identical to Equations (4) and (22). The numerical scheme described by Equation (23) was presented by Rees et al. (1989) under the name of DELO and subsequently re-baptized by De la Cruz Rodríguez & Piskunov (2013) as DELO-linear. It is a second-order accurate method, as shown in Figure 2. The explicit values of the coefficients ${{\boldsymbol{\Phi }}}_{k}$, ${{\boldsymbol{\Phi }}}_{k+1}$, ${{\boldsymbol{\Psi }}}_{k}$, and ${{\boldsymbol{\Psi }}}_{k+1}$ are provided in Appendix B.

The natural successive step is a parabolic Lagrangian interpolation of the effective source function (Murphy 1990), i.e., Equation (20) with q = 2. Considering three spatial points $\{{\tau }_{k-1},{\tau }_{k},{\tau }_{k+1}\}$ located along the optical depth grid, the effective source function ${\bf{S}}$ is approximated by the parabolic interpolation,

$$\begin{eqnarray*}{\bf{S}}(\tau ,{\boldsymbol{I}}(\tau )) & \approx & {{\bf{S}}}_{k+1}\displaystyle \frac{(\tau -{\tau }_{k})(\tau -{\tau }_{k-1})}{{\rm{\Delta }}{\tau }_{k}({\rm{\Delta }}{\tau }_{k}+{\rm{\Delta }}{\tau }_{k-1})}\\ & & -\ {{\bf{S}}}_{k}\displaystyle \frac{(\tau -{\tau }_{k+1})(\tau -{\tau }_{k-1})}{{\rm{\Delta }}{\tau }_{k}{\rm{\Delta }}{\tau }_{k-1}}\\ & & +\ {{\bf{S}}}_{k-1}\displaystyle \frac{(\tau -{\tau }_{k+1})(\tau -{\tau }_{k})}{{\rm{\Delta }}{\tau }_{k-1}({\rm{\Delta }}{\tau }_{k}+{\rm{\Delta }}{\tau }_{k-1})},\end{eqnarray*}$$

for τ in the interval $[{\tau }_{k},{\tau }_{k+1}]$. The necessity of a third interpolation point cannot be satisfied by numerical values at ${\tau }_{k+2}$, because no information is available for ${{\boldsymbol{I}}}_{k+2}$. After two integrations by parts of Equation (18) and some algebra, one obtains

$$\begin{eqnarray}&&{{\boldsymbol{\Phi }}}_{k+1}{{\boldsymbol{I}}}_{k+1}={{\boldsymbol{\Phi }}}_{k}{{\boldsymbol{I}}}_{k}+{{\boldsymbol{\Phi }}}_{k-1}{{\boldsymbol{I}}}_{k-1}+{{\boldsymbol{\Psi }}}_{k+1}+{{\boldsymbol{\Psi }}}_{k}+{{\boldsymbol{\Psi }}}_{k-1}.\end{eqnarray} \tag{ 24 }$$

The method described by Equation (24), which might be called DELO-parabolic, is third-order accurate, as shown in Figure 2. The coefficients ${{\boldsymbol{\Phi }}}_{k-1}$, ${{\boldsymbol{\Phi }}}_{k}$, ${{\boldsymbol{\Phi }}}_{k+1}$, ${{\boldsymbol{\Psi }}}_{k-1}$, ${{\boldsymbol{\Psi }}}_{k}$, and ${{\boldsymbol{\Psi }}}_{k+1}$ are provided in Appendix B. The higher order of accuracy is essential to detect, for instance, the second-order behavior of the emission vector often present in realistic atmospheric models, avoiding its possible systematic overestimation.

DELO-constant and DELO-linear formal solvers compute the Stokes ${{\boldsymbol{I}}}_{k+1}$ solely on the basis of information about the preceding Stokes value ${{\boldsymbol{I}}}_{k}$ and are classified as one-step methods. In this sense, one-step methods have no memory, i.e., they forget all of the prior information that has been gained. In contrast, DELO-parabolic takes into account the two most recently found Stokes vectors ${{\boldsymbol{I}}}_{k}$ and ${{\boldsymbol{I}}}_{k-1}$, entering in the class of the multistep methods.

This family of formal solvers can be further expanded by just increasing the interpolation degree q of the effective source function. For instance, a third-order Lagrangian polynomial would generate DELO-cubic. However, the complexity of the numerical methods would increase as well and the adaptation to non-uniform grids would become gradually more cumbersome.

One should emphasize that the detailed behavior of the error curves plotted in Figures 2 and 3 depends on the considered atmospheric model. Indeed, the main scope of these figures is to highlight the overall order of accuracy of the various methods. It would certainly be wrong to try to reach conclusions on the performance of different methods on the basis of a qualitative comparison of small details of the error curves, obtained considering a single model atmosphere. The atmospheric model considered in this work is described in Appendix C.

Zoom In Zoom Out Reset image size

The choice of the Lagrangian form for the polynomial ${{\boldsymbol{P}}}_{q}$ in Equation (17) is certainly not univocal and the literature provides different interpolation methods. Mihalas et al. (1978) used the Hermitian interpolation for the integration of the scalar radiative transfer equation. An interesting set of suitable interpolants was proposed by Auer (2003) for the same problem: among them the monotonic Hermite interpolants recently used by Ibgui et al. (2013) in the IRIS code. In the same year, De la Cruz Rodríguez & Piskunov (2013) applied Bézier polynomials to the DELO strategy, generating the quadratic and cubic DELO-Bézier methods.

A detailed description of these methods, and a comparison with other methods, such as DELO-linear and DELOPAR, can be found in the above-mentioned publications, and will not be repeated here. As shown in Figure 3, both quadratic and cubic DELO-Bézier methods show fourth-order accuracy. It is interesting to observe that when treating smooth functions, the Bézier curves introduced by De la Cruz Rodríguez & Piskunov (2013) are forced to be identical to the Hermite polynomials of corresponding degree by adopting very specific control points. A detailed discussion of Hermitian methods will be presented in the second paper of this series, which will focus on high-order methods.

Aiming to increase the order of accuracy with respect to the DELO-linear strategy, Trujillo Bueno (2003) opted for a semi-parabolic interpolation, namely a parabolic interpolation for the modified emission vector $\tilde{{\boldsymbol{\epsilon }}}$ and a linear interpolation for the ${\bf{K}}{\boldsymbol{I}}$ term. The parabolic interpolation of $\tilde{{\boldsymbol{\epsilon }}}$ is performed by considering the three spatial points $\{{\tau }_{k},{\tau }_{k+1},{\tau }_{k+2}\}$, which differs from the set used by DELO-parabolic. In the case of a diagonal dominant propagation matrix, this strategy can increase the order of accuracy in the Stokes parameter I, provided a high-order integration of the opacity (see Appendix A). However, the local truncation error in the Stokes Q, U, and V, is still dominated by the linear approximation of the term ${\bf{K}}{\boldsymbol{I}}$, resulting in a second-order method as shown in Figure 3. Therefore, the lack of improvement with respect to DELO-linear comes directly from the design of the method and not, as erroneously conjectured, from the so-called overshooting. This should also explain the unsatisfactory performance of DELOPAR found by De la Cruz Rodríguez & Piskunov (2013), as presented in their error curves. The DELOPAR strategy remains a honest method for the non-polarized case or for vanishing dichroism and anomalous dispersion coefficients. In fact, in both cases, the modified propagation matrix ${\bf{K}}$ disappears and a parabolic approximation of the source term $\tilde{{\boldsymbol{\epsilon }}}$, supported by a proper conversion to optical depth (see Appendix A), produces an effective third-order numerical scheme. This can be appreciated in Štěpán & Trujillo Bueno (2013), where the log–log error figure clearly shows the superiority of DELOPAR over DELO-linear for the scalar case. Štěpán & Trujillo Bueno (2013) applied a similar technique to create BESSER, the formal solver used in the PORTA code. The same argumentation can be applied to it, explaining its second-order accuracy confirmed by Figure 3.