Research Note
Broyden’s method for the solution of the multilevel nonLTE radiation transfer problem
Institut de Recherche en Astrophysique et Planétologie, Université de
Toulouse, CNRS,
14 av. E. Belin,
31400
Toulouse,
France
email: fpaletou@ast.obsmip.fr
Received:
13
October
2010
Accepted:
12
December
2010
This study concerns the fast and accurate solution of multilevel nonLTE radiation transfer problems. We propose and evaluate an alternative iterative scheme to the classical MALI method. Our study is instead based on the application of Broyden’s method for the solution of nonlinear systems of equations. Comparative tests, in 1D planeparallel geometry, of the popular MALI method and our alternative method are discussed. The Broyden method is typically 4.5 times faster than MALI. This makes it also fairly competitive with the GaussSeidel and Successive OverRelaxation methods developed after MALI.
Key words: radiative transfer / methods: numerical
© ESO, 2011
1. Introduction
The solution of the nonLTE multilevelatom radiative transfer problem is a classical one in astrophysics. Indeed, consistently with the departure of the source functions from Planck functions, the assumption of nonLTE 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 nonLTE case, one has on the contrary to solve simultaneously and selfconsistently 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 deexcitation 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 Λiterationbased 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 SocasNavarro & Trujillo Bueno (1997).
In this study, we investigate the use of a quasiNewton 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 WolfRayet stars. Broyden’s method was more recently invoked in the context of the coupledescape 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 speedup, 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 deexcitation rates. In general, these rates depend on the electronic density, so that if the latter is not known a priori, terms like n_{i}C_{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 wellknow shortcharacteristics method with monotonic parabolic interpolation introduced by Auer & Paletou (1994).
In that frame, and using the ShermanMorrison 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_{0}, at every τ depth, and an initial Broyden matrix B_{0}; then we compute . The iterative scheme is (7)then update (8)then compute (9)and finally, update (10)while ∥F∥ _{2} > ε. In practice, we chose ε = 10^{2}, which guarantees that we reached a fully converged state (see also Fig. 1).
Fig.1 Typical relative error, R_{k} (thin), convergence error, C_{k} (thick) and F_{k} (dotted) for MALI (dash) and Broyden (full) vs. the number of iterations for both schemes. 

Open with DEXTER 
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 modelatom. In this way, the grand function F defined by Eqs. (5) and (6) can be fully evaluated. Then, we compute an initial Jacobian, B_{0}, 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_{0}.
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 semiinfinite plane parallel slab model of τ_{max} = 10^{10}, discretised in a various number of points per decade in optical depth, using also a 3 to 10 energy level Hatom, 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^{5} s^{1}. We also adopted the definitions initially proposed by Auer et al. (1994) for the relative error, from one iteration (k) to another:
and for the “convergence error”:
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
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, C_{k}, for MALI (dashed) and Broyden (full) are displayed vs. the computing time for a 5level H atom and five (thin) and height (thick) points per decade in optical depth. 

Open with DEXTER 
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 5level 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 wellconverged 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_{0}. 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 Hatom model, for MALI (dashed) and Broyden (full). The dashdotted curve corresponds to the time required for the evaluation of the initial Jacobian. 

Open with DEXTER 
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 nonLTE radiative transfer problem were introducted by the GaussSeidel (GS) and successive overrelaxation (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 selfconsistent evaluation of the electron density in a multilevel nonLTE problem, if necessary – a problem for which MALI additionally needs a NewtonRaphson 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 nonLTE 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 wellknown MALI method. Another advantage is that it does not require any modification of usual formal solvers, as is the case for GSSOR methods developed after MALI. It is also potentially very wellsuited for parallel computing. Further tests will include the selfconsistent treatment of the ionization balance, which is usually treated together with MALI with a NewtonRaphson scheme. In a next step, we will consider more demanding models such as H_{2}O, for instance.
References
 Auer, L. H. 1991, in Stellar Atmospheres: Beyond Classical Models, NATO ASI Series (Dordrecht: Reidel) [Google Scholar]
 Auer, L. H., & Mihalas, D. 1969, ApJ, 158, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Auer, L. H., & Paletou, F. 1994, A&A, 285, 675 [NASA ADS] [Google Scholar]
 Auer, L. H., Fabiani Bendicho, P., & Trujillo Bueno, J. 1994, A&A, 292, 599 [NASA ADS] [Google Scholar]
 Broyden, C. G. 1965, Math. Comp., 19, 577 [CrossRef] [Google Scholar]
 Elitzur, M., & Asensio Ramos, A. 2006, MNRAS, 365, 779 [NASA ADS] [CrossRef] [Google Scholar]
 Heinzel, P. 1995, A&A, 299, 563 [NASA ADS] [Google Scholar]
 Koesterke, L., Hamman, W.R., & Kosmol, P. 1992, A&A, 255, 490 [NASA ADS] [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres (San Francisco: Freeman) [Google Scholar]
 Paletou, F. 1995, A&A, 302, 587 [NASA ADS] [Google Scholar]
 Paletou, F., & Léger, L. 2007, JQSRT, 103, 57 [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., et al. 1992, Numerical recipes (Cambridge: University Press) [Google Scholar]
 Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171 [NASA ADS] [Google Scholar]
 SocasNavarro, H., & Trujillo Bueno, J. 1997, ApJ, 490, 383 [NASA ADS] [CrossRef] [Google Scholar]
 Trujillo Bueno, J., & Fabiani Bendicho, P. 1995, ApJ, 455, 646 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig.1 Typical relative error, R_{k} (thin), convergence error, C_{k} (thick) and F_{k} (dotted) for MALI (dash) and Broyden (full) vs. the number of iterations for both schemes. 

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

Open with DEXTER  
In the text 
Fig.3 Respective computing times vs. number of levels of the Hatom model, for MALI (dashed) and Broyden (full). The dashdotted curve corresponds to the time required for the evaluation of the initial Jacobian. 

Open with DEXTER  
In the text 