© ESO, 2011

1. Introduction

The solution of the non-LTE multilevel-atom radiative transfer problem is a classical one in astrophysics. Indeed, consistently with the departure of the source functions from Planck functions, the assumption of non-LTE implies that the population density of the atomic or molecular levels considered departs from what can be derived at LTE in a straightforward manner with Saha and Boltzmann relations (see e.g., Mihalas 1978).

In the non-LTE case, one has on the contrary to solve simultaneously and self-consistently for a set of N_T equations of radiative transfer together with N_L equations of statistical equilibrium (hereafter ESE) that describe the detailed balanced of excitation and de-excitation processes between every atomic or molecular levels. Because absorption and stimulated emission radiative rates explicitely depend on the radiation field, which itself depends on the level populations, this problem is intrinsically a search for the solution of coupled nonlinear equations.

Since the beginning of numerical radiative transfer in the late 60’s, the two most popular methods used for tackling this problem have been the complete linearization method of Auer & Mihalas (1969) and the accelerated Λ-iteration-based scheme called MALI (Rybicki & Hummer 1991). Despite their apparent differences, they however have in common that basically they deal with linearized equations. An interesting comparative study of these two approaches has been made by Socas-Navarro & Trujillo Bueno (1997).

In this study, we investigate the use of a quasi-Newton numerical method for the solution of the nonlinear ESE. Our choice is Broyden’s method (1965), whose elements will be presented in Sect. 2. To the best of our knowledge, Koesterke et al. (1992) were the first to bring this numerical scheme into the field of radiation transfer. Their study was presented in the context of the modelling of spherically expanding atmospheres of hot and massive Wolf-Rayet stars. Broyden’s method was more recently invoked in the context of the coupled-escape probability method (Elitzur & Asensio Ramos 2006).

Besides the required algebra and the mention of caveats with regard to the implementation of the method, it remains difficult however to figure out from Koesterke et al. (1992) the actual performance of this approach. The authors also briefly compared their method with another one, but they mentioned a significant speed-up, which is provided by Broyden algorithm for large N_L atomic models. In particular, it is a pity that no comparison with the MALI method was made yet, because it is contemporary with the publication of Rybicki & Hummer (1991). This evaluation is the aim of the present work.

2. The numerical scheme

For an N_L-level atomic model, the ESE will are generally written as a set of elementary equations: (1)where the A_ij and B_ij stand for the spontaneous emission respectively, and the absorption and stimulated emission rates, n_i represents the population density for each energy level, and is the scattering integral for each radiatively allowed transition we will consider. Besides the radiative processes, the C_ij are collisional excitation and de-excitation rates. In general, these rates depend on the electronic density, so that if the latter is not known a priori, terms like n_iC_ij are nonlinear in the population densities. Below we will consider only cases for which the collisional rates are known a priori.

The scattering integral entering the ESE is formally written as (2)where, assuming complete redistribution in frequency, the source function is defined as (3)A large system of Eqs. (1) is homogeneous unless one of the equations is replaced by a constraint equation like, for instance, a conservation equation of the form (4)

2.1. Broyden’s algorithm

The system of Eqs. (1) and (4) can be reformulated by defining a function F that acts on the set of , where τ is a discrete depth along the opacity scale used to sample our slab or atmosphere. F is defined as (5)for i ≠ N_L and, if i = N_L(6)Computations of F require repeated evaluations of the scattering integrals – Eq. (2) – which we perform here using the well-know short-characteristics method with monotonic parabolic interpolation introduced by Auer & Paletou (1994).

In that frame, and using the Sherman-Morrison formula (see e.g., Press et al. 1992), which provides an analytical formula for the direct computation of the inverse of the Broyden matrix, our algorithm consists of the following steps. We choose an initial vector n₀, at every τ depth, and an initial Broyden matrix B₀; then we compute . The iterative scheme is (7)then update (8)then compute (9)and finally, update (10)while  ∥F∥ ₂ > ε. In practice, we chose ε = 10^-2, which guarantees that we reached a fully converged state (see also Fig. 1).

Fig.1 Typical relative error, Rk (thin), convergence error, Ck (thick) and Fk (dotted) for MALI (dash) and Broyden (full) vs. the number of iterations for both schemes.

2.2. Initialization

The proper initialization of the Broyden scheme is a critical issue. We employed the following method, which was tested as suitable both from the standpoint of an adequate start and from the one of an acceptable computing time.

Before starting the iterative scheme, we assume LTE populations for our model-atom. In this way, the grand function F defined by Eqs. (5) and (6) can be fully evaluated. Then, we compute an initial Jacobian, B₀, using the finite difference scheme fdjac described by Press et al. (1992).

In the figure below with regard to the timing properties of Broyden’s method vs. MALI, we will always include the specific time necessary for the evaluation of B₀.

3. Comparison of Broyden vs. MALI

We adopted the popular “flavour” of the MALI method with a diagonal approximate operator as described by Rybicki & Hummer (1991) and without acceleration of convergence schemes.

In order to compare the properties of the Broyden scheme with those of MALI, we adopted a 1D semi-infinite plane parallel slab model of τ_max = 10¹⁰, discretised in a various number of points per decade in optical depth, using also a 3 to 10 energy level H-atom, inspired by the classic benchmark proposed by Avrett (1968; see also Léger & Paletou 2007). As in the latter, the slab temperature was fixed at 5000 K and the collisional rates set at 10⁵   s^-1. We also adopted the definitions initially proposed by Auer et al. (1994) for the relative error, from one iteration (k) to another:

(11)

and for the “convergence error”:

(12)

where n_(∞) is the “fully converged” solution obtained for a given method and model after a substantial number of iterations. We also introduce the quantity

(13)

i.e., the Euclidian norm of F, the function defined by Eqs. (5) and (6). Note that for the MALI method is defined by a modified Eq. (5) following the preconditioning strategy proposed by Rybicky & Hummer (1991).

Fig.2 Convergence error, Ck, for MALI (dashed) and Broyden (full) are displayed vs. the computing time for a 5-level H atom and five (thin) and height (thick) points per decade in optical depth.

3.1. Convergence

In Fig. 1 we display the rates of convergence of the Broyden and MALI methods. The convergence error C_k for each method was computed with population densities obtained once F_k < 10^-2 in both cases. The case used here is a 5-level H atom with a 1D slab discretized by a 5 points per decade grid in optical depth. It is worth noting that to reach F_k < 10^-2 and a well-converged solution, one should iterate until one reaches R_k values as small as 10^-10 typically. In terms of the number of iterations, Broyden typically beats MALI by more than one order of magnitude. However, the quite distinct nature of each method makes this comparison incomplete. Below we perform this analysis but we will compare the respective computing times.

3.2. Sensitivity to the spatial (optical depth) refinement

In Fig. 2 we turn to an analysis of the respective timing properties of Broyden and MALI. We show that Broyden, again, always beats MALI by a typical factor of the order of about 4–5 in time. This is also the case when the spatial grid refinement is increased from 5 to 8 points per decade, for instance. It is important to note that timings given for Broyden include the computation of the initial matrix B₀. This is why the rates of convergence displayed for the Broyden method do not start at t = 0 in Fig. 2.

3.3. Sensitivity to the number of transitions

A next important point to investigate regards the advantage of Broyden against MALI when an increasing number of atomic transitions is considered. Again, as demonstrated in Fig. 3, the full Broyden iterative process is always significantly faster than MALI. In general, the gain in total computing time with the Broyden method is of the order of 4–5. This is less than the the gain of the order of 8 already reported by Koesterke et al. (1992), although their method of reference was presumably different from MALI.

Fig.3 Respective computing times vs. number of levels of the H-atom model, for MALI (dashed) and Broyden (full). The dash-dotted curve corresponds to the time required for the evaluation of the initial Jacobian.

3.4. Discussion

We are aware that the MALI method can be accelerated by specific schemes (see e.g., Auer 1991). But the most significant improvements in the field of iterative methods for the non-LTE radiative transfer problem were introducted by the Gauss-Seidel (GS) and successive over-relaxation (SOR) methods (Trujillo Bueno & Fabiani Bendicho 1995). It was already shown, for instance, that SOR always beats both Jacobi (i.e., accelerated Λ-iteration with the diagonal of the full Λ operator as an approximate operator) and GS, even when Ng acceleration of convergence is applied.

Beyond the fact that Broyden is significantly faster than MALI, we can also add that Broyden is as competitive as the SOR method, according to Paletou & Léger (2007; see their Table 1 where comparable timing and the corresponding iteration numbers were given for MALI, as we used in the present study, GS and SOR).

The Broyden method is also potentially more advantageous than MALI and GS/SOR because of its intrinsic capability to deal with the self-consistent evaluation of the electron density in a multilevel non-LTE problem, if necessary – a problem for which MALI additionally needs a Newton-Raphson scheme, as proposed by Heinzel (1995) and Paletou (1995).

Another important point is that, as indicated in our Fig. 3, a great deal of time of our Broyden code is spent on the computation of the initial Jacobian, a task which can be performed with great advantage using parallel computing. The inner structure of the fdjac routine indeed permits parallelization with a high scalability.

As a final comment, it is also important to consider that Broyden’s method can be easily implemented to already existing codes without the need of modifying the formal solver, unlike the GS/SOR methods.

4. Conclusion

We propose an alternative method for the solution of the non-LTE multilevel radiation transfer problem. It is based on Broyden’s method for the solution of nonlinear systems of equations. The method is easy to implement and is about of factor of 4.5 times faster than the well-known MALI method. Another advantage is that it does not require any modification of usual formal solvers, as is the case for GS-SOR methods developed after MALI. It is also potentially very well-suited for parallel computing. Further tests will include the self-consistent treatment of the ionization balance, which is usually treated together with MALI with a Newton-Raphson scheme. In a next step, we will consider more demanding models such as H₂O, for instance.

References

All Figures

	Fig.1 Typical relative error, Rk (thin), convergence error, Ck (thick) and Fk (dotted) for MALI (dash) and Broyden (full) vs. the number of iterations for both schemes.
In the text

	Fig.2 Convergence error, Ck, for MALI (dashed) and Broyden (full) are displayed vs. the computing time for a 5-level H atom and five (thin) and height (thick) points per decade in optical depth.
In the text

	Fig.3 Respective computing times vs. number of levels of the H-atom model, for MALI (dashed) and Broyden (full). The dash-dotted curve corresponds to the time required for the evaluation of the initial Jacobian.
In the text