1 Introduction

The radiative transfer equation for polarized light (RTE) has resisted all the approaches for a practical general solution even in those cases in which the source function is well determined (Landi Degl'Innocenti & Landi Degl'Innocenti 1985). In the familiar scalar case when polarization is not taken into account, the solution of the RTE is straightforward, provided that the source function is given as a function of the integration variable (i.e. optical or geometrical depth). This is not the case, however, for polarized light. Going from one to four dimensions introduces new and fundamental problems. Paper I demonstrated that the non-commutativity properties of a general absorption matrix are at the origin of all those difficulties, and proposed a way of dealing with them. The insights gained from this work have indeed suggested a method to obtain a general analytical solution not based on infinite series (Magnus 1954; Landi Degl'Innocenti 1987). This solution has been presented by López Ariste & Semel (1999) but further work will be necessary before a numerical code based on it can be implemented. On the other hand, the solution found in Paper I for the case of commuting absorption matrices proved to be very powerful, and at the same time very easy to implement into a numerical code. The set of commuting matrices includes all the cases for which a particular solution is already known (Unno-Rachkowsky solution - Rachkowsky 1967-, constant matrix case -Landi Degl'Innocenti & Landi Degl'Innocenti 1985-, ...) as well as some new non-trivial cases, for instance those in which the azimuth of magnetic field varies linearly with depth. In order to extend this approach to more general absorption matrices, including non-commuting ones, it was proposed in Paper I to firstly approximate any absorption matrix by a commuting one, for which a solution was known, and secondly to improve this solution by linearizing the residual. In this paper we will present the numerical code that implements that algorithm, which we have called DIAGONAL, and will study its performances using DELO (Rees et al. 1989) as a benchmark.

When introducing a new numerical code, analysis of errors and convergence rates is very important. In relation to RTE, this analysis has once again been hindered by the absence of an adequate analytical solution. Comparison with Rachkowsky solution (Rachkowsky 1967) in a Milne-Eddington atmosphere with constant absorption matrix, used to serve as a first indication of how good an algorithm was or how it converged (Rees et al. 1989; Bellot Rubio et al. 1998). This is however not a conclusive test because, due to its linear source function, any integration method more accurate than first order in the integration variable (usually z the geometrical light path, or $\tau$ the optical depth) will give an exact solution for these kinds of atmosphere models. Although other solutions with non-linear source functions (Lites et al. 1988) or non-constant absorption matrices (Landi Degl'Innocenti & Landi Degl'Innocenti 1985) are known, they have apparently never been used for this purpose. Therefore performances have been usually calculated by comparison to previous algorithms, and so we shall do the same with DIAGONAL. The basic assumption under this comparison test is that any accepted algorithm gives the exact solution (up to a given numerical precision) for a large number of integration layers. This assumption is revealed to be perfectly well founded. Nevertheless, the study of the algorithms when compared to an analytical solution for non-trivial cases seems to be very desirable. The algorithm proposed in Paper I provides an exact solution for some special non-trivial cases and, once the numerical code is shown to perform correctly, this fact can be used to test other numerical algorithms.

In the same line of thought we will make use in Sect. 2 of the formal solution given by Magnus (1954, see also the Appendix of Paper I for a resumée of this seminal paper) to characterize the algorithm in terms of convergence and errors. In Sect. 3 we will show in detail the way in which the proposed algorithm has been implemented. Finally Sect. 4 will be dedicated to show some numerical comparisons of the new code, called DIAGONAL, with the well-known DELO, and, as proposed above, of the convergence of DELO to the correct solution given by DIAGONAL in particular (but not trivial) cases.