# Formal solutions for polarized radiative transfer I. the DELO family

## Abstract

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

###### Subject headings:

Radiative transfer – Polarization – Methods: numerical## 1. Introduction

The transfer of partially polarized light is described by a system of coupled first-order, inhomogeneous ordinary differential equations. Explicitly,

(1) |

where is the spatial coordinate measured along the ray under consideration, is the Stokes vector,

is the propagation matrix, and is the emission vector, which represents the source term. The different coefficients appearing in the propagation matrix and in the emission vector depend on the considered frequency, on the propagation direction, and on different atmospheric parameters (Landi Degl’Innocenti & Landi Degl’Innocenti, 1985; Landi Degl’Innocenti, 1987). For notational simplicity, the frequency dependence of the quantities is not explicitly indicated.

Analytical solutions of Equation (1) are available for a few simple model atmospheres only (López Ariste & Semel, 1999a), which explains the necessity of numerical schemes able to solve it. The definition of *formal solution* was first introduced for the scalar problem: it is the evaluation of the radiation intensity, given knowledge of the boundary conditions and the spatial, angular, and frequency dependence of the opacity and the emissivity at a discrete set of points (Mihalas, 1978; Auer, 2003). The generalization to the polarized case consists in substituting radiation intensity, opacity, and emissivity by Stokes vector, propagation matrix, and emission vector, respectively.

The formal solution of the radiative transfer equation is a key step of iterative schemes for solving the nonlinear full radiative transfer problem, where the atomic system and the radiation field interact in non-local thermodynamic equilibrium (non-LTE) conditions. From the computational point of view, the formal solution is, in most cases, the slowest part of the iterative scheme (e.g., Štěpán & Trujillo Bueno, 2013). Moreover, large magnetohydrodynamic simulations of stellar atmospheres call for massive synthesis of Stokes profiles. Therefore, the requirement for the numerical method is to be as accurate and as fast as possible. The effort of the community has produced an extensive literature on the different formal solvers, the major contributions being summarized in Table 1. The quest of the “best formal solver” available is still open and the comparison between the different numerical schemes provided by the community is somehow confusing.

This paper aims to give a structured overview over formal solvers, in particular over those belonging to the Diagonal Element Lambda Operator (DELO) family, and to clarify some incoherences found in the literature. Section 2 presents a reference paradigm for the characterization of formal solvers. The concepts of *order of accuracy*, *stability* and *computational cost* are briefly presented and used to characterize a reference numerical scheme, namely the trapezoidal method. Section 3 provides an introduction to exponential integrators, a class of numerical methods for the solution of differential equations. A simple description of this class is given, in an attempt to highlight its specific features and its paternity of the DELO methods. In Section 4, the well known DELO family is presented and characterized. Some new methods belonging to this family are introduced and compared to the already existing ones. Finally, Section 5 provides remarks and conclusions, with a view on future work.

Year |
Method |
Proposed by |
---|---|---|

1974 | Runge-Kutta-Merson | Wittmann (1974) |

1976 | Runge-Kutta 4 | Landi Degl’Innocenti (1976) |

1985 | Piecemeal Evolution Operator | Landi Degl’Innocenti & Landi Degl’Innocenti (1985) |

1989 | Zeeman Feautrier | Rees, Durrant, & Murphy (1989) |

1989 | DELO-linear | Rees, Durrant, & Murphy (1989) |

1998 | (cubic) Hermitian | Bellot Rubio, Ruiz Cobo, & Collados (1998) |

1999 | DIAGONAL | López Ariste & Semel (1999b) |

2003 | DELOPAR | Trujillo Bueno (2003) |

2013 | (quadratic and cubic) DELO-Bézier | De la Cruz Rodríguez & Piskunov (2013) |

2013 | BESSER | Štěpán & Trujillo Bueno (2013) |

2016 | Piecewise Continuous | Steiner, Züger, & Belluzzi (2016) |

## 2. Characterization of formal solvers

In order to be able to answer the question “which is the best formal solver?”, one first has to fix some criteria for judging the different numerical schemes. Briefly discussing a method, claiming some strong points, and showing them with a few specific examples, is not the best way to proceed. The literature about the numerical approach to ordinary differential equations is particularly broad and great efforts have been directed toward the characterization of the different methods. This section aims to give an overview on this characterization, in order to facilitate the comprehension of the advantages, weaknesses, and possible incoherences of the already existing formal solvers and of those yet to come. To ease the appreciation, a very common and simple method is presented and analyzed, namely the trapezoidal method.

### 2.1. Exempli gratia: the trapezoidal method

A spatial grid is introduced, discretizing the ray path. The spatial coordinate and the index increase along the propagation direction. For a given grid point , the points and represent the *upwind* and *downwind* points, respectively, guaranteeing . Applying the fundamental theorem of calculus to Equation (1) in the interval , one obtains

(2) |

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

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

(3) |

where is the cell width. The numerical approximation of a certain quantity at node is indicated by substitution of the explicit dependence on with the subscript , for instance,

Inserting numerical quantities for , , and in Equation (3) and applying some algebra, one recovers the implicit linear system

(4) |

where the matrices and and the vectors and are given by

Given the upwind Stokes vector at node , one solves the implicit linear system given by Equation (4), finding the emergent Stokes vector at node .

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.

### 2.2. Order of accuracy

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

supplied by the discrete grid , the local truncation error is defined by

where and are the exact solutions at nodes and , respectively, and represents the numerical approximation at node after performing the single step from to . The operator represents any suitable vector norm. A numerical method for an ordinary differential equation has order of accuracy if the local truncation error, which usually depends on the step size denoted by , satisfies

(5) |

Hence, the larger the order of accuracy, the faster the error is reduced as decreases.

By contrast, the *global error* is defined as

and represents the accumulation of the local truncation error over all the 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 implies a global error of . Explicitly

(6) |

using Equation (5) and . The constant typically depends on the exact solution and on other parameters of the numerical scheme, but is strictly independent of . 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 must be replaced by the maximal step size given by

with . The expression “suitable discrete grid” indicates that the maximal step size is inversely proportional to the total number of grid points, i.e., .

The power law presented in Equation (6) implies a linear relationship between the logarithms of the global error and the step size , i.e.,

where . 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 . The trapezoidal method is well known to be second-order accurate (e.g., Frank & Leimkuhler, 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 in order to avoid the *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 is dominated by the first term for small enough , but for large step sizes higher-order terms may tangibly contribute, depending on their constants and .

### 2.3. Stability

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

(7) | ||||

with . The solution converges to zero as for . A Runge-Kutta method applied to the IVP (7) can be recast into the form

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 and this condition is equivalent to

(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 and the eigenvalue , more precisely to the term . The *stability region* of a Runge-Kutta method is defined as the set of complex values 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

(9) | ||||

where the matrix has a basis of eigenvectors corresponding to the eigenvalues . 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 & Leimkuhler (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 . It is important to mention the fact that the eigenvalues of the propagation operator 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

The stability region for the trapezoidal method is then given by the condition (8) and it is usually displayed as presented in Figure 1a. If the stability region of a numerical method contains the whole left-hand side of the complex plane, as in the case of Figure 1a, 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 (Collatz, 1966; Dahlquist, 1963), but a deep digression would stray from the main aim of this paper.

A strong limitation of this simplified stability analysis is the assumption of a constant eigenvalue in Equation (7). A less restrictive 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 reads

(10) |

where and are the eigenvalues at the positions and , respectively. Thus, the stability region depends on both eigenvalues and , 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 , supporting the assumption of a constant eigenvalue in the stability analysis.

### 2.4. Computational cost

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 directly from the precedent value , i.e.,

while implicit one-step methods find by solving an equation of the type

which results in the additional solution of a 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.

## 3. Exponential integrators

Exponential integrators form a class of numerical methods for ordinary differential equations. This class is based on the exact integration of the linear part of the IVP, aiming at reducing the stiffness of the differential equation.

In order to present this class, one considers the following general IVP

(11) | ||||

which is equivalent to Equation (1) given the initial value . One splits into linear and nonlinear contributions, i.e.,

where the matrix does not depend on the variable and the nonlinear term is given by . Equation (11) is then recast in the form

(12) | ||||

and the exact integration of the linear part in the interval yields the *variation of constants* formula

(13) |

The integral in Equation (13) has to be numerically approximated and a large variety of different options is available, e.g., Runge-Kutta discretizations as explained by Cox & Matthews (2002). Moreover, for a non-diagonal matrix , the evaluation of the matrix exponential often requires an approximation and one has to combine the integrator with well-chosen algorithms from numerical linear algebra.

In the following section, a specific strategy is applied, where only the nonlinear term is approximated and the exponential operator is treated exactly.

## 4. The DELO family

In this section, a particular family of methods belonging to the class of exponential integrators is presented. The first method of this family applied to the formal solution for polarized light was proposed by Rees et al. (1989) under the name of DELO. Thereafter, a second version by Trujillo Bueno (2003) took the appellative DELOPAR. Additional improvements, in terms of Bézier interpolations, were recently provided by De la Cruz Rodríguez & Piskunov (2013) and by Štěpán & Trujillo Bueno (2013).

As exhaustively explained by Guderley & Hsu (1972), this technique takes into account analytically the diagonal elements of the propagation matrix , aiming to remove stiffness from the problem. Therefore, the well-known radiative transfer equation for polarized light given by Equation (1) is brought in the form given by Equation (12), with the additional constraint of a diagonal matrix . This reformulation is facilitated by the fact that the diagonal elements of the propagation matrix are all identical. Replacing the coordinate by the optical depth defined by

(14) |

one recasts Equation (1) into

(15) |

where represents the identity matrix. The quantity is the *effective source function* and it is defined by

with the modified propagation matrix (whose diagonal elements are all equal to zero) and the modified emission vector .

Observe that the optical depth scale is in fact defined so that the coordinate , defined in Equation (14), decreases along the ray path, thus . Providing the upwind Stokes , the operator on the left of Equation (15) is inverted leading to the formula

(16) |

which is analogous to Equation (13).

As anticipated, different numerical quadratures of the integral in Equation (16) lead to different numerical schemes. In particular, the DELO technique approximates the effective source function by a polynomial of degree inside the integration interval, i.e.,

(17) |

Observing that and evaluating Equation (16) in the interval , one obtains

(18) |

where . The integral can then be solved by parts, yielding an implicit or explicit linear system for the Stokes vector .

As explained by Guderley & Hsu (1972), the local truncation error is due to the fact that the effective source function is approximated by a polynomial of degree . According to Henrici (1962), this approximation results in the following local truncation error

(19) |

and, from definition (5), a DELO method involving a polynomial should show an order of accuracy equal to . Equation (19) is not defined for , i.e., for a constant approximation of the effective source function. In this case, the local truncation error corresponds to the one for , following a behavior similar to the implicit midpoint method (Deuflhard & Bornemann, 2002). Moreover, the evaluation of the quantity plays a fundamental role in the accuracy of the numerical scheme, as explained in Appendix A.

As already mentioned, the DELO strategy is thought to remove stiffness from the problem. In fact, a method based on Equation (16) tends to A-stability for a vanishing modified propagation matrix. Therefore, a sufficiently small matrix should imply a rather wide stability limit. However, the simple stability analysis presented in the previous section cannot be applied to the present family of methods, because of the two different contributions from the exponential terms and the modified propagation matrix. Nevertheless, some indications can be deduced if the simpler scalar case is analyzed. Guderley & Hsu (1972) show that the DELO strategy increases stability with respect to the corresponding Adams-Bashford methods (Deuflhard & Bornemann, 2002).

On the other hand, the second-order accurate trapezoidal method already guarantees A-stability, which dispenses from a stability improvement for formal solvers having an order of accuracy .

### 4.1. DELO-constant, DELO-linear and DELO-parabolic

In a very general way, the polynomial in Equation (17) is obtained by the Lagrangian interpolation

(20) |

where and the Lagrange basis polynomials are given by

Here 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

(21) |

where the different coefficients and depend on the chosen polynomial and on the numerical values and . Provided the previous Stokes for , the linear system (21) can be solved to obtain the numerical approximation at . Therefore, once given the boundary condition at , the recursive application of Equation (21) provides the emergent Stokes vector at the end of the ray path.

As first choice, the effective source function is assumed constant inside the interval , i.e.,

and can be approximated by the midpoint rule

Replacing the polynomial 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.,

(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 , , , and are provided in Appendix B.

The next step is to obtain a relation between and when approximating the effective source function by a linear interpolation, i.e., Equation (20) with ,

for in the interval . One proceeds with an analytical integration and after some algebra obtains

(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 , , , and 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 . Considering three spatial points located along the optical depth grid, the effective source function is approximated by the parabolic interpolation,

for in the interval . The necessity of a third interpolation point cannot be satisfied by numerical values at , because no information is available for . After two integrations by parts of Equation (18) and some algebra, one obtains

(24) |

The method described by Equation (24), which might be called DELO-parabolic, is third-order accurate, as shown in Figure 2. The coefficients , , , , , and 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 solely on the basis of information about the preceding Stokes value 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 and , entering in the class of the multistep methods.

This family of formal solvers can be further expanded by just increasing the interpolation degree 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.

### 4.2. DELO-Bézier methods

The choice of the Lagrangian form for the polynomial 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.

### 4.3. Particular cases: DELOPAR and BESSER

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 and a linear interpolation for the term. The parabolic interpolation of is performed by considering the three spatial points , 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 , provided a high-order integration of the opacity (see Appendix A). However, the local truncation error in the Stokes , , and , is still dominated by the linear approximation of the term , 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 disappears and a parabolic approximation of the source term , 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.

Formal solver |
Order of accuracy |
---|---|

DELO-constant | 2 |

DELO-linear | 2 |

DELO-parabolic | 3 |

DELOPAR | 2 |

BESSER | 2 |

Quadratic DELO-Bézier | 4 |

Cubic DELO-Bézier | 4 |

## 5. Conclusions

This paper pays particular attention to the characterization of the different formal solvers for polarized radiative transfer. The paradigmatic analysis for numerical schemes proposed here is based on three different criteria: *order of accuracy*, *stability*, and *computational cost*.

The order of accuracy of a numerical scheme indicates how fast the error decreases, when reducing the step size. Therefore, it can be only appreciated by considering the global (or local) error dependence on the step size and not through single Stokes profiles.

When discussing numerical methods, the term stability relates to the attenuation of numerical errors with the recursive application of a scheme and not to the erratic behavior of high-order polynomial approximation (the so-called overshooting). One must always guarantee that any possible instability is avoided, when judging the order of accuracy of a method.

The computational cost is related to the complexity of the algorithm. Therefore, one suggests to maintain as bony a method as possible, avoiding unnecessary superstructures.

In this perspective, some considerations about the DELO family can be exposed. Regarding the order of accuracy, one realizes that DELO methods do not converge better than more conventional methods: DELO-constant, DELO-linear, and DELOPAR are only second-order accurate, the two-step DELO-parabolic method effectively reaches third-order accuracy, and quadratic and cubic DELO-Bézier methods usually perform as fourth-order accurate methods, as summarized in Table 2. One must point out that second-order accuracy is already guaranteed by the simple trapezoidal method. The discussion about stability is somehow more involved. The DELO approach is thought to remove stiffness from the problem, but no stability improvement is necessary for formal solvers having an order of accuracy , because the trapezoidal method already guarantees A-stability as shown in Figure 1a. However, one could still promote the DELO strategy for high-order methods (). Concerning the computational cost, DELO methods differ from standard implicit numerical methods for two reasons: the DELO strategy requires an additional conversion to optical depth (see Appendix A), which, however, is anyway important to mitigate the variation of the eigenvalues of the propagation operator along the ray path, enforcing stability. Moreover, DELO coefficients require the evaluation of exponential terms (see Appendix B) and this results in extra computational cost, as explained in Section 2.4.

In conclusion, the necessity of the DELO strategy for the numerical treatment of the polarized radiative transfer is not warranted for low-order methods and not well motivated for high-order methods. A second paper, focused on high-order formal solvers, will try to build a clear hierarchy with respect to order of accuracy, stability and computational cost. The effective performances when dealing with realistic atmospheric models remain to be explored.

## Appendix A Conversion to optical depth

As already pointed out by De la Cruz Rodríguez & Piskunov (2013), many radiative transfer applications require a conversion of the spatial scale, e.g., from geometrical height to optical depth . From Equation (14) one obtains

(A1) |

The numerical integration introduces an error, which could lead to a reduced order of accuracy of the formal solver. In practice, a trapezoidal integration of Equation (A1) is inadequate to perform numerical schemes based on high-order interpolations of the effective source function (e.g., DELO-parabolic). Therefore, high-order DELO methods require a corresponding high-order numerical evaluation of the integral in Equation (A1).

## Appendix B DELO coefficients

In order to keep the notation as close as possible to Rees et al. (1989), the following definitions are introduced

The coefficients of the DELO-constant method, Equation (22), are given by

The coefficients of the DELO-linear method, Equation (23), are given by

The coefficients of the DELO-parabolic method, Equation (24), are given by

with

The DELO-parabolic coefficients could suffer of problematic division with vanishingly small quantities. Thus, in case of small , a Taylor expansion of the exponential term to third-order is indicated.

## Appendix C Atmospheric model

The atmosphere model and parametric description used for the calculations shown in this paper are very similar to the ones used by Steiner et al. (2016). The radiative transfer is computed at different frequencies in a spectral interval containing a hypothetical, magnetically sensitive spectral line. The problem is formulated in reduced frequencies

(C1) |

with the line-center frequency and the Doppler width. The spectral interval is considered. The atmosphere is assumed to be plane-parallel, and the chosen reference spatial coordinate is the continuum optical depth along the vertical direction, defined by

where is the continuum absorption coefficient at the line-center frequency, and is the geometrical height (increasing in the outward direction). The atmosphere extends in the range .

No scattering or atomic polarization is taken into account, so that polarization is only introduced by the Zeeman effect. The magnetic field vector and the absorption and anomalous dispersion profiles (which enter in the definition of the coefficients of the propagation matrix and of the emission vector ) are assumed to be depth-independent. Under these assumptions, the propagation matrix can be parametrized in the form (see, e.g., Landi Degl’Innocenti & Landolfi, 2004; Steiner et al., 2016)

where is the ratio between the frequency-integrated line absorption coefficient and the continuum absorption coefficient, and is a constant matrix. The emission vector, on the other hand, can be written as

where is the usual intensity source function. The following analytical form of and (the only depth-dependent quantities of the problem) have been considered

(C2) | ||||

(C3) |

The exponential term in Equation (C2) is included in order to reproduce a possible emission rise at low optical depths (e.g., at chromospheric heights). The results of Figures 2 and 3 have been obtained on the basis of the smooth variation of and shown in Figure 4, with the following values of the parameters appearing in Equations (C2) and (C3): , , , , , and .

The intensity of the magnetic field is specified through the dimensionless parameter , being the Larmor frequency. The results shown in Figures 2 and 3 have been obtained for , and assuming the magnetic field to have an inclination with respect to the vertical. A damping constant has been chosen for calculating the Voigt and Faraday-Voigt profiles entering the definition of the radiative transfer coefficients. No macroscopic (bulk) velocity has been considered. The calculations have been performed for the radiation propagating outwards in the atmosphere, along the vertical direction. At the bottom of the atmosphere (boundary condition) an unpolarized radiation field has been introduced. The Stokes parameters of the emergent radiation have been calculated in units of the parameter , whose exact value is thus irrelevant for the calculations shown in this paper. The reference direction for positive Stokes has been taken in the plane defined by the vertical and by the magnetic field.

## Appendix D Error calculation

Denoting with and the reference and the numerically computed emergent Stokes vectors, respectively, at the reduced frequency , the global error for the th Stokes vector component is computed as

(D1) |

The error is given by the maximal discrepancy between the reference and the simulated Stokes parameter over the entire profile, normalized by the maximal amplitude in the reference profile. Equation (D1) is not defined for a constant profile, because of a vanishing denominator. In that case, one needs to introduce a different error definition. The reference emergent Stokes profile is the exact solution approximated by means of high-order numerical methods, using a hyperfine grid sampling with more than points-per-decade of the continuum optical depth. Different high-order methods (e.g., DELO-parabolic and quadratic DELO-Bézier) are used to cross-check the reference emergent profile.

### Footnotes

- affiliationmark:
- affiliationmark:
- affiliationmark:
- affiliationmark:
- affiliationmark:
- affiliationmark:
- affiliationmark:

### References

- Auer, L. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 288, Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, 3
- Bellot Rubio, L. R., Ruiz Cobo, B., & Collados, M. 1998, ApJ, 506, 805
- Collatz, L. 1966, The numerical treatment of differential equations (Springer-Verlag)
- Cox, S., & Matthews, P. 2002, J. Comput. Phys., 176, 430
- Dahlquist, G. G. 1963, BIT Numerical Mathematics, 3, 27
- De la Cruz Rodríguez, J., & Piskunov, N. 2013, ApJ, 764, 33
- Deuflhard, P., & Bornemann, F. 2002, Scientific Computing with Ordinary Differential Equations (Springer-Verlag)
- Frank, J., & Leimkuhler, B. 2012, Computational Modelling and Dynamical Systems, Lecture Notes (Universities of Amsterdam and Edinburgh)
- Goldreich, O. 2008, Computational Complexity: A Conceptual Perspective, 1st edn. (Cambridge University Press)
- Guderley, G., & Hsu, C.-C. 1972, Math. Comp., 26, 51
- Hackbusch, W. 2014, The Concept of Stability in Numerical Mathematics (Springer Berlin Heidelberg)
- Henrici, P. 1962, Discrete variable methods in ordinary differential equations (Wiley)
- Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013, A&A, 549, A126
- Landi Degl’Innocenti, E. 1976, A&A, 25, 379
- Landi Degl’Innocenti, E. 1987, in Numerical Radiative Transfer, ed. W. Kalkofen (Cambridge University Press), 265–278
- Landi Degl’Innocenti, E., & Landi Degl’Innocenti, M. 1985, Sol. Phys., 97, 239
- Landi Degl’Innocenti, E., & Landolfi, M. 2004, Astrophysics and Space Science Library, Vol. 307, Polarization in Spectral Lines (Kluwer Academic Publishers)
- López Ariste, A., & Semel, M. 1999a, A&A, 350, 1089
- —. 1999b, A&AS, 139, 417
- Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (W.H. Freeman and Company)
- Mihalas, D., Auer, L. H., & Mihalas, B. R. 1978, ApJ, 220, 1001
- Murphy, G. A. 1990, PhD thesis, , Univ. Sidney, (1990)
- Rees, D. E., Durrant, C. J., & Murphy, G. A. 1989, ApJ, 339, 1093
- Schörghofer, N. 2015, Lessons in Scientific Computing (University of Hawaii at Manoa)
- Steiner, O., Züger, F., & Belluzzi, L. 2016, A&A, 586, A42
- Trujillo Bueno, J. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 288, Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, 551
- Štěpán, J., & Trujillo Bueno, J. 2013, A&A, 557, A143
- Wittmann, A. 1974, Sol. Phys., 35, 11