2 Error analysis

Usually (e.g. Paper I, López Ariste & Semel 1999 or the original paper by Landi Degl'Innocenti & Landi Degl'Innocenti 1985), the formal solution to the RTE is presented as  

$$\vec{I}(z)=\int_{z_0}^z \tens{O}(z,z')\vec{J}(z'){\rm d}z' + \tens{O}(z,z_0)\vec{I}(z_0) ,$$ (1)

where $\vec{I}$ is the Stokes vector, $\vec{J}$ the emission vector, z the light path through the atmosphere and $\tens{O}(z,z_0)$ the so-called evolution operator. This expression can be seen as the particular solution of an inhomogeneous first-order differential equation whose homogeneous counterpart has a solution given by the evolution operator:

$$\frac{{\rm d}}{{\rm d}t}\tens{O}(t,z_0)=-\tens{K}(t) \tens{O}(t,z_0),$$ (2)

where $\tens{K}$ denotes the absorption matrix and where

$$\tens{O}(z_0,z_0)={\mathchoice {\rm 1\mskip-4mu l} {\rm 1\mskip-4mu l} {\rm 1\mskip-4.5mu l} {\rm 1\mskip-5mu l}},$$

the $4\times 4$ identity matrix. Magnus (1954) provides us with a general solution to this kind of equation in the form of a unique exponential:

$$\tens{O}(z,z_0)={\mbox e}^{-\Omega (z,z_0)}.$$

The exponent $\Omega (z,z_0)$ is given as an infinite series whose 2 first terms are  

$$\Omega (z,z_0)=\int _{z_0}^z \tens{K} (t) {\rm d}t + \frac{1... ...t _{z_0}^t \tens{K} (\sigma) d\sigma \right] {\rm d}t + \ldots$$ (3)

The rest of the terms in the series are higher recursive commutators between the absorption matrix and its integral. This solution, impractical from a numerical point of view, will allow us nevertheless to describe and clarify the basic sources of numerical errors when integrating the RTE. Introducing Eq. (3) in Eq. (1) we obtain:  

$$\vec{I}(z)=\int_{z_0}^z {\mbox e}^{-\Omega (z,z')}\vec{J}(z'){\rm d}z' + {\mbox e}^{-\Omega (z,z_0)}\vec{I}(z_0).$$ (4)

Any numerical integration algorithm must solve, implicitly or explicitly, three problems:

1.

Determination of the exponent $\Omega (z,z')$.

2.

Calculation of the exponential.

3.

Integration of the emission vector.

Below we shall discuss these three issues in greater detail.

2.1 Determination of the exponent $\Omega (z,z')$

The only known attempt to explicitly calculate the second or higher terms in Eq. (3) was done by López Ariste & Semel (1998) for a matrix linear with the integration variable. It was given for illustrative purposes and reflects the difficulties of this computation in the form given by Magnus. The analytical solution given by López Ariste & Semel (1999) gives some hope as to the inclusion of the commutators of the absorption matrix with its integral in the solution, but some further work is still necessary. From the numerical codes side, it is usually difficult to disentangle the determination of $\Omega (z,z')$ from the problem of calculation of the exponential. For instance, DELO integrates the diagonal of the absorption matrix and linearizes all the rest. This linear approximation should take into account the calculation of the exponential of the off-diagonal terms as well as the higher terms of $\Omega (z,z')$. In practice, DELO imposes a constant matrix for each layer, the variations with depth of the model atmosphere being taken into account by the computation of numerous layers. For a constant matrix, all the commutators are null and $\Omega (z,z')$ becomes very easy to calculate. The Canarias Hermite method (Bellot Rubio et al. 1998; Ruiz Cobo et al. 1998) uses a completely different reasoning: instead of worrying about the formal solution of the RTE, this method never integrates explicitly the equation but takes into account the derivatives of the solution, which can be obtained directly from the differential equation. In that sense this method implicitly calculates the whole Magnus' series, as it is built such as to ensure that its derivatives satisfy always the RTE. On the other hand, the number of derivatives which can be calculated and the number of points at which they can be calculated is limited. The result is equivalent to a sampling of the solution, with the usual associated errors of any sampling. Perhaps the most important among them are the errors associated with the calculation of the exponential (see next subsection) and the problem of the oscillating evolution operator (Bellot Rubio et al. 1998), which will be briefly discussed below.

The solution proposed in Paper I, and implemented in DIAGONAL, integrates explicitly the RTE. The absorption matrix is separated into two pieces: one ($\tens{K}_{\rm o}$) which commutes with its own integral and another which does not commute. For the first one $\Omega (z,z')$ is immediately determined from the first term of the Magnus' series:

$$\Omega_{\rm o} (z,z') = \int_z^{z'} \tens{K}_{\rm o} (t) {\rm d}t.$$

The difference with the constant-matrix solution is that Paper I provides the necessary and sufficient conditions for any absorption matrix to commute. DIAGONAL takes advantage of that and can, for instance, incorporate some gradients into the commuting part. The best example is the one of a linear variation with optical depth of the azimuth of magnetic field (see Paper I).

  

Figure 1: Stokes profiles for the 6301 Å  Fe I line in a model with magnetic field equal to $2300 + 800 \log{\tau}$ G, inclination equal to $45^{\circ}+10\log{\tau}$, velocity in the l.o.s. equal to $0.3\log{\tau}$ km s-1, constant azimuth equal to $40^\circ$ and ratio of the line to continuum absorption coefficient equal to 40

For the residual part we made a choice and we linearize it in a similar way to DELO with the off-diagonal terms. The residual is very small compared to the whole, so that this linearization is a fair approximation. Furthermore, the atmosphere is divided therefore improving the linear approximation and diminishing the error. It is important to remark that the number of layers necessary to keep under control errors is very small compared to DELO, as most of the problem is solved analytically through $\Omega _{\rm o}$.

  

Figure 2: Linear correction in percent of the total signal applied to I and V of Fig. 1 with 15 integration layers

The choice of a DELO-like linear approximation is not the sole one available and even not the best. The residual part obeys Eq. (66) of Paper I and a variety of numerical methods can be thought of to solve it. Our choice responds to simplicity. As we shall see later, the transformations involved in the analytical part of the solution put a price in terms of computation time per layer. Further sophistications in solving the residual risked increasing this time and diminishing the performance of the code.

How good is the hypothesis that the residual is really a residual, i.e., small compared with the analytical part? To answer this question we have drawn the residual in terms of percentage of the final signal at each wavelength calculated and at each node in depth, for a model with strong gradients in intensity and azimuth of magnetic field and in velocity in the l.o.s. (see the next section for a full characterization of the model atmosphere used in this and other tests presented in this paper). The Stokes profiles produced by this model atmosphere for the Fei 6301 Å  line are presented in Fig. 1. Figures 2 and 3 show the residual when the model is integrated using 15 layers.

  

Figure 3: Linear correction in percent of the total signal applied to Q and U of Fig. 1 with 15 integration layers

In view of them, both the name residual applied to the non-commuting part of the initial absorption matrix and the linear approximation applied to it are completely justified. Of course this residual correction becomes smaller and smaller as we approach the surface: the layers are distributed logarithmically and they are thinner near $\tau = 0$. The thinner the layer the better is the approximation of a commuting matrix for it. It is also interesting to note the presence of ripples in the correction surface for Stokes I and V. They are due to the oscillations in the evolution operator. These oscillations are one of the hardest problems for any integrator based on a sampling of the solution. If the sampling is not adequated to the frequency of oscillation an error appears. This problem has been discussed at length by Bellot Rubio et al. (1998). Any analytical solution overrides this problem, but the residual represented in Figs. 2 and 3 is based on a 2-point sampling (the linear approximation), and hence the problem arises.

The number of layers used in the previously referred figures is higher than the one usually used with DIAGONAL. Fewer layers mean thicker layers so that the correction becomes larger. We present in Figs. 4 and 5 the residuals for the same model atmosphere of the previous example integrated this time with just 3 layers. Although the corrections are bigger, the linear approximation continues to be completely satisfactory even in these extreme cases.

  

Figure 4: Linear correction in percent of the total signal applied to I and V of Fig. 1 with 3 integration layers

  

Figure 5: Linear correction in percent on the total signal applied to Q and U of Fig. 1 with 3 integration layers

2.2 Calculation of the exponential

Surprisingly most algorithms fail to solve this calculation, quite anodyne in the scalar case. Once more the extrapolation from the scalar case is not evident. The exponential of a matrix is defined as  

$${\mbox e}^{\Omega} = {\mathchoice {\rm 1\mskip-4mu l} {\rm 1... ...+ \Omega + \frac{1}{2}\Omega ^2+\frac{1}{3!} \Omega ^3 + \ldots$$ (5)

and it is in general not easy to calculate the entire series. Although convergence can be ensured, the series is not a power expansion. Only the presence of the growing factorial in the denominator provides a hint on how many terms may suffice to attain a given precision, and in any case the calculation is too slow. The different numerical methods used for the integration of the RTE stop the series at different points: Runge-Kutta method (Landi Degl'Innocenti 1976) stops at the 4th order, DELO at the 2nd order (Rees et al. 1989) and the Canarias Hermite method attain the 4th order too. In view of the convergence rates of all these methods, it shall be not too far from the truth to say that the calculation of the exponential is the first source of errors in all those methods (instability problems aside).

In fact, at the cost of some more algebra, other efficient methods can be used to calculate the exponential of a matrix. We will mention here two, perhaps the most evident. The first one is to calculate the exponential of the associated diagonal matrix: if $\lambda _i$ are the eigenvalues of $\Omega$, and $\tens{T}$ the diagonalization matrix, then

$${\mbox e}^{\Omega}= \tens{T} \pmatrix{{\mbox e}^{\lambda _1}... ...lambda _3} & 0\cr0&0&0 &{\mbox e}^{\lambda _4}} \tens{T}^{-1}.$$

This calculation requires knowledge of the eigenvalues and of the eigenvectors (to calculate $\tens{T}$). It is a very interesting method whenever we deal with an analytical matrix, whose expression is known in advance. Analytical values for the eigenvalues and eigenvectors can be calculated once for all, and the numerical code just works out the particular case. This is the actual method used in DIAGONAL: the formal absorption matrix is well-known and indeed its eigenvalues and eigenvectors have already been calculated in Paper I. Hence, the exponential is calculated exactly by DIAGONAL at the price of one matrix multiplication (this interesting particularity of DIAGONAL is the origin of its name).

A second method goes through the Hamilton-Cayley theorem particularized to the exponential function. It is known from linear algebra that any matrix satisfies its characteristic polynomial. For a $4\times 4$ matrix this is a 4th order polynomial equation:

$$a \Omega ^4 + b\Omega ^3 +c \Omega ^2 + d\Omega + e = 0.$$

We can write $\Omega ^4$ as a function of the other powers of the matrix, and do it for any power higher than 4. The exponential series in (5) becomes a polynomial of order 3. The correct values of the coefficients of this polynomial are given by the Hamilton-Cayley theorem. This method is particularly advised for numerical matrices for which the diagonal method cannot be efficiently employed.

2.3 Integration of the emission vector

The integration of the inhomogeneous part of the transfer equation can be addressed by two main methods. The first one is to continue with the sampling technique as the Canarias Hermite method does. We have discussed it in the subsection on the calculation of the exponent. The second one is to give an analytical expression for the emission vector. For instance, we can assume it is linear with optical depth, as in the Milne-Eddington model atmosphere. Once we have an analytical function we can integrate the whole analytically or use numerical integration techniques. In the case of DIAGONAL we have chosen the analytical integration and used for all the tests a Milne-Eddington linear source function. This is a choice made for those tests. Other functions can be used without changing the characteristics of the code, although it is worthy to note that the analytical integration becomes rapidly very difficult due to the presence of an exponential of non-linear functions. In any case, errors do not accumulate due to this integration in the actually used code.

2.4 Summary of errors

As a conclusion of this section, we can summarize the sources of errors in the DIAGONAL code by saying that at worst it will behave like DELO, i.e., like a 2nd order integration technique. DIAGONAL cannot be qualified as an nth order technique in general. The reason is its mixed character. If it is asked to integrate a commuting absorption matrix it will provide an analytical result. In the worst case, when the residual becomes too big, DIAGONAL will behave as DELO. Nevertheless, this last extreme has never been found in the actual tests.